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Abstract 

This  study  investigates  the  use  of  a  linear  quadratic  terminal  controller  to  recon¬ 
figure  satellite  formations  using  atmospheric  drag  actuated  control  while  minimizing  the 
loss  of  energy  of  the  formation.  The  linearized  Clohessy-Wiltshire  equations  of  motion  are 
used  to  describe  the  motion  of  the  two-satellite  formation  about  an  empty  reference  posi¬ 
tion  maintained  at  the  formation  center.  Reconfigurations  to  final  in-plane  and  elliptical 
formations  are  simulated  at  orbital  radii  of  6800  km  and  7000  km,  and  the  altitude  loss 
and  a  Av  budget  were  recorded  as  performance  measures  for  each  reconfiguration.  The 
final  states  of  the  spacecraft  upon  reconfiguration  were  propagated  forward  in  time  over 
20  orbital  periods  to  ensure  the  final  conditions  were  achieved.  Simulations  proved  that 
minimizing  the  loss  of  orbital  energy  effectively  minimizes  the  loss  in  altitude,  and  drag 
actuated  control  is  fully  capable  of  controlling  the  radial  and  in-track  motion  of  satellite 
formations,  although  the  cross-track  motion  is  uncontrollable. 
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SATELLITE  FORMATION  CONTROL 


USING  ATMOSPHERIC  DRAG 

I.  Introduction 

1.1  Satellite  Formation  Overview 

As  technology  evolves,  new  ideas  develop  to  improve  current  processes.  This  evolu¬ 
tion  is  found  in  all  systems,  including  the  space  arena.  Satellite  formation  flying  is  one 
area  of  improvement  for  satellite  systems  that  has  gained  considerable  interest  in  recent 
years.  Although  the  idea  has  been  around  for  over  two  decades,  the  implementation  of 
formation  flying  is  still  in  its  infancy.  As  a  result,  research  continues  to  develop  in  hopes 
of  improving  the  level  of  confidence  and  understanding  of  satellite  formation  flying. 

A  satellite  formation  consists  of  two  or  more  satellites  flying  in  close  proximity  to 
one  another  in  a  specified  geometry,  cooperating  as  a  distributed  system  to  achieve  a  given 
mission  [6].  This  distributed  nature  of  satellite  formations  provides  several  advantages 
over  a  single  satellite  system,  leading  to  potential  increases  in  capability  for  space  systems. 
One  area  of  increased  capability  presented  by  formation  flying  is  remote  sensing.  Remote 
sensing  systems,  such  as  radar,  are  limited  by  the  aperture  size  of  the  antenna  on  the 
satellite.  This  aperture  is,  in  turn,  hindered  by  physical  limitations  since  the  antenna  must 
either  be  small  enough  to  launch  into  orbit  whole  or  launched  in  pieces  and  constructed 
on-orbit,  a  difficult  endeavor  to  say  the  least.  Moreover,  the  control  and  stability  of 
large  optical  systems  further  restrict  the  current  capabilities.  Satellite  formations,  on  the 
other  hand,  offer  greater  angular  and  spatial  resolution  through  the  creation  of  large, 
synthetic  apertures,  formed  by  precisely  positioned  satellites  equipped  with  distributed 
optical  systems  [6].  The  aperture  could  then  be  adjusted  as  desired  by  increasing  or 
decreasing  the  size  of  the  formation  or  by  reconfiguring  the  formation,  alleviating  the 
physical  restrictions  on  the  aperture. 

Aside  from  increased  capabilities  in  remote  sensing,  satellite  formations  provide  po¬ 
tential  cost  savings,  increased  reliability,  and  flexibility.  The  cost  savings  are  directly  tied 
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to  the  complexity  of  the  system.  Due  to  the  importance  and  substantial  financial  invest¬ 
ment  involved  with  satellite  systems,  spacecraft  are  often  pushed  to  the  physical  size  limits 
prescribed  by  the  launch  vehicles  as  a  result  of  redundant  systems  placed  onboard  to  help 
reduce  the  risk  of  mission  failure.  By  replacing  large,  complex,  single  satellite  systems  with 
a  cluster  of  small,  simply  designed  satellites  carrying  distributed  payloads  and  cooperating 
to  complete  the  mission,  the  system  can  yield  substantial  savings  in  cost  and  liability.  For 
example,  should  one  satellite  in  the  cluster  fail,  the  other  satellites  in  the  formation  can 
still  be  used  to  at  least  partially,  if  not  completely,  accomplish  the  mission.  The  formation 
can  be  easily  reconfigured  to  optimize  its  capability,  and,  if  needed,  additional  satellites 
can  be  launched  to  restore  the  full  capability  of  the  system.  Furthermore,  as  technology 
develops  or  the  mission  changes,  satellite  formations  can  easily  be  upgraded  with  replace¬ 
ment  satellites  as  the  mission  sees  fit.  This  redundancy  and  flexibility  provides  a  potential 
advantage  in  the  military  arena  as  well.  With  the  development  of  anti-satellite  technology, 
a  distributed  system  could  prove  more  resistant  to  attack  since  an  adversary  would  have 
to  destroy  multiple  satellites  to  cripple  the  fleet,  not  just  one  [5].  The  system  flexibility 
is  seen  in  the  development  of  the  system  on  the  ground  as  well.  Very  often,  the  launch  of 
satellites  are  delayed  due  to  setbacks  during  procurement  or  development,  many  times  due 
to  a  single  component.  For  a  distributed  system,  the  individual  satellites  can  be  launched 
as  they  are  developed,  allowing  partial  performance  of  the  system  until  all  of  the  satellites 
are  on-orbit,  thereby  mitigating  any  critical  delays. 

Although  satellite  formations  provide  extensive  advantages  over  a  single  satellite  sys¬ 
tem,  key  issues  that  are  not  as  prevalent  for  a  single  satellite  can  quickly  complicate  the 
system.  One  key  issue  is  maintaining  the  desired  formation,  or  formation  keeping.  As 
discussed  by  Sabol  et  ah,  several  factors  can  lead  to  formation  degradation  [18].  Semi¬ 
major  axis  differences  result  in  different  orbital  periods,  causing  the  formation  to  diverge. 
Eccentricity  and  inclination  differences  will  cause  the  orbital  planes  to  separate  as  an  effect 
of  J'2 ,  a  topic  that  will  be  discussed  further  in  Section  2.1.4.  Furthermore,  differences  in 
mean  anomaly  can  lead  to  an  in-track  drift.  These  perturbing  factors  are  all  related  to 
orbital  element  differences,  but  physical  differences  between  the  satellites  themselves  may 
result  in  formation  degradation  as  well.  Orbital  perturbations  from  Earth’s  full  geopo- 
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tential,  third-body  effects,  atmospheric  drag,  solar  radiation  pressure,  and  other  natural 
effects  can  cause  a  formation  to  separate  if  the  satellites  are  not  identical.  As  an  example, 
a  satellite  with  a  larger  presented  area  per  mass  will  likely  experience  more  atmospheric 
drag,  resulting  in  a  higher  in-track  deceleration,  causing  the  formation  to  drift  apart.  Con¬ 
sequently,  having  an  idea  of  these  effects  is  essential  to  the  design  of  the  formation  and  to 
the  formation  scheme  selected. 

Two  common  formation  designs  are  the  leader-follower  formation  and  the  empty 
reference  position  formation.  The  leader-follower  formation  consists  of  a  leader  satellite  in 
a  reference  orbit  and  one  or  more  follower  satellites  whose  position  and  motion  are  defined 
relative  to  the  leader.  This  formation  is  ideal  since  it  provides  a  physical  reference  point  for 
the  follower  satellites.  However,  one  potential  drawback  from  the  design  is  that  the  follower 
satellites  are  usually  following  a  different  orbit  than  the  leader,  causing  the  formation  to 
naturally  degrade.  The  empty  reference  position  formation,  on  the  other  hand,  consists  of 
two  or  more  satellites  positioned  with  respect  to  some  reference  position  at  the  center  of 
the  formation,  unoccupied  by  another  satellite.  This  design  lends  itself  well  to  formation 
flying  since  all  of  the  satellites  in  the  formation  usually  have  very  similar  orbits,  reducing 
the  degradation  of  the  formation  from  orbital  perturbations.  These  two  formation  designs 
will  be  further  elaborated  upon  in  Section  2.1.4. 

Regardless  of  the  formation  scheme  selected,  each  satellite  formation  must  be  capable 
of  performing  a  reconfiguration,  whether  the  maneuver  is  required  for  formation  keeping 
purposes  or  to  transfer  into  a  new  formation  to  perform  some  aspect  of  the  mission.  The 
most  common  and  practical  means  of  formation  control  is  through  the  use  of  conventional 
thrusters.  The  technology  involving  chemical  thrusters  is  well  understood,  and  their  high 
impulse  nature  allow  researchers  to  treat  orbital  maneuvers  as  velocity  change  impulses, 
simplifying  the  dynamics  of  the  transfer.  However,  one  drawback  to  conventional  thrusters 
is  the  requirement  of  onboard  fuel.  Onboard  expendables  increase  the  size  and  weight  of 
the  satellite  and  necessitates  strict  fuel  conservation  to  prolong  the  life  of  the  formation. 

Another  technique  for  formation  control  is  a  passive  control  approach  using  elec- 
trodynamic  tethers.  This  technique  involves  a  formation  of  satellites  that  are  attached 
to  one  another  via  electrodynamic  tethers.  A  current  is  then  passed  through  the  tether 
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to  control  the  formation.  Since  the  current  generated  through  the  tether  is  produced  by 
electricity  from  fuel  cells  or  solar  power,  the  expendables  required  onboard  the  spacecraft 
are  minimal,  reducing  the  overall  cost  of  development  and  launch,  and  possibly  increasing 
the  lifetime  of  the  formation  since  the  effective  life  is  no  longer  dictated  by  the  amount  of 
fuel  onboard.  But  tethers  also  present  obvious  difficulties  that  must  be  addressed.  Since 
the  satellites  are  connected,  tangling  the  tethers  during  reconfiguration  must  be  avoided. 
Also,  should  a  tether  sever  due  to  high  stress,  fatigue,  or  errant  space  debris,  the  means 
for  control  is  lost. 

A  third  technique  that  has  been  considered  for  formation  control  is  atmospheric 
drag.  Using  atmospheric  drag  for  control  is  desirable  because,  like  the  electrodynamic 
tethers,  drag  does  not  require  fuel  onboard  the  spacecraft,  and  it  is  in  infinite  supply  for 
all  near-Earth  satellites.  By  simply  rotating  either  the  solar  panels  or  drag  panels  attached 
to  the  satellite,  the  amount  of  drag  acting  on  the  satellite  can  be  varied,  providing  a 
source  of  control.  However,  exposing  a  satellite  to  increased  drag  also  causes  the  orbit  to 
degrade  more  rapidly.  Without  adequate  fuel  onboard  the  satellite  to  boost  the  orbit,  this 
accelerated  rate  of  decay  will  shorten  the  vehicle’s  lifetime,  requiring  altitude  conservation 
methods  to  extend  the  life  of  the  formation. 

1.2  Literature  Review 

The  control  of  satellite  formations  using  conventional  thrusters  has  been  extensively 
studied.  Ulybyshev  used  a  discrete-time  linear  quadratic  regulator  (LQR)  feedback  con¬ 
trol  law  to  perform  station  keeping  of  a  low-Earth  orbit  satellite  constellation  perturbed 
by  atmospheric  drag  [21].  The  control  law  minimized  both  the  along-track  relative  dis¬ 
placements  and  the  orbital  period  differences  between  the  satellites  to  maintain  the  relative 
position  of  the  satellites,  proving  a  linear-quadratic  control  law  can  effectively  handle  the 
control  of  multiple  satellites  for  formation  keeping.  Irvin  compared  linear  and  nonlinear 
control  techniques  for  the  reconfiguration  of  a  satellite  formation  [5].  LQR,  LQR  with  lin¬ 
earizing  feedback,  and  the  state  dependent  Riccati  equation  (SDRE)  control  methods  were 
applied  to  both  the  linear  and  nonlinear  forms  of  the  Clohessy-Wiltshire  equations,  and  the 
performance  of  each  control  method  was  compared  to  a  near-optimal,  open-loop,  discrete- 
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time,  impulsive  maneuver.  The  study  found  that  the  linear  controllers  performed  better 
than  the  nonlinear  control  law,  but  the  continuous  feedback  controller  was  less  optimal 
than  the  discrete-time,  impulsive  maneuver  control.  Sparks  compared  the  formation  keep¬ 
ing  control  costs  for  a  LQR  feedback  control  law  formulated  from  the  linear,  unperturbed 
equations  of  motion  to  the  control  costs  to  directly  overcome  the  gravitational  perturba¬ 
tions  acting  on  a  formation  [19] .  The  LQR  control  law  was  simulated  using  a  high-fidelity 
model  incorporating  gravitational  perturbations,  atmospheric  drag,  and  other  perturbing 
effects.  The  study  found  that  as  the  time  interval  between  burns  is  increased,  the  error 
between  the  actual  and  desired  formations  increase,  but  the  control  costs  approach  the 
theoretical  minimum  to  overcome  the  gravitational  perturbations.  Sabol  et  al.  looked  at 
formation  keeping  costs  of  four  basic  formations  of  two  satellites:  in-plane,  in-track,  cir¬ 
cular,  and  projected  circular  [18].  The  in-plane  formation  consists  of  the  two  satellites 
flying  in  the  same  orbital  plane,  but  separated  by  mean  anomaly.  The  in-track  formation 
is  designed  such  that  the  satellites  share  the  same  ground  track,  placing  the  satellites  in 
different  orbital  planes  separated  by  right  ascension  of  the  ascending  node  to  account  for 
the  Earth’s  rotation.  The  circular  formation  is  designed  such  that  the  satellites  maintain  a 
set  distance  from  one  another,  making  the  follower  appear  to  orbit  the  reference  satellite. 
And  the  projected  circular  formation  is  designed  such  that  the  satellite  cluster  maintains 
a  fixed  distance  in  the  in-track/cross-track  plane.  The  study  used  a  high  fidelity  perturba¬ 
tion  model  to  propagate  the  orbital  elements  of  the  satellites  forward  in  time  to  determine 
the  stability  of  the  formation  and  the  formation  keeping  control  costs.  All  results  yielded 
rather  high  control  costs  to  achieve  the  desired  conditions. 

While  extensive  research  has  been  performed  in  the  area  of  formation  keeping  using 
conventional  thrusters,  very  little  work  has  been  done  using  atmospheric  drag  for  formation 
control.  Leonard  looked  at  using  differential  drag  between  two  spacecraft  for  formation 
keeping  using  a  two-part  control  law  [10].  The  main  control  law  drove  the  average  position 
of  the  satellite  to  a  target  position  relative  to  the  reference  satellite.  A  separate,  second 
control  law  was  then  used  to  minimize  the  eccentricity  of  the  orbit  while  maintaining  that 
average  position.  Both  satellites  were  equipped  with  drag  panels  that  could  be  oriented 
at  either  0  or  90  degrees  to  provide  a  set  positive,  negative,  or  zero  relative  acceleration. 
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This  technique  yielded  favorable  results  as  the  satellite  was  driven  to  the  desired  position 
in  less  than  24  hours  for  all  cases  tested,  and  the  eccentricity-minimizing  control  law  was 
able  to  maintain  the  position  within  4.17  feet. 

More  recently,  Wedekind  considered  formation  control  using  differential  drag  on  a 
leader- follower  formation  [23].  In  this  case,  a  proportional-plus-integral  controller  was  used 
to  continuously  alter  the  cross-sectional  area  of  the  follower  satellite  to  counter  a  drag  bias 
between  the  two  satellites  in  the  formation.  The  control  of  three  different  formations, 
in-plane,  in-track,  and  circular,  was  considered.  Wedekind  achieved  favorable  results  for 
these  three  formations  when  the  drag  biases  were  kept  small,  showing  that  a  formation 
drag  bias  can  be  mitigated  and  the  formation  maintained  by  altering  the  presented  area 
of  the  follower  satellite. 

1.3  Spacecraft  Formation  Missions 

While  extensive  research  into  formation  flying  is  being  performed,  several  systems 
are  being  developed  and  launched  to  prove  the  viability  of  such  techniques.  NASA  has 
embraced  the  idea  of  smaller  and  cheaper  spacecraft,  filling  numerous  missions  to  validate 
technology  and  prove  formation  flying  concepts  [6].  NASA’s  Earth  Observing-1  (EO-1)  mis¬ 
sion,  launched  in  2001,  “developed  and  validated  a  number  of  instrument  and  spacecraft 
bus  breakthrough  technologies  designed  to  enable  the  development  of  future  earth  imaging 
observatories  that  will  have  a  significant  increase  in  performance  while  also  having  reduced 
cost  and  mass”  [14].  EO-1  was  inserted  into  a  formation  with  the  Landsat  7  satellite  to 
validate  autonomous  operations  concepts  and  other  technologies  onboard  the  spacecraft 
as  part  of  NASA’s  “smaller,  cheaper  and  more  capable”  spacecraft  tests.  NASA’s  Space 
Technology  5  mission,  launched  in  2006,  demonstrated  the  enhanced  capabilities  forma¬ 
tion  flying  provides  to  scientific  research  while  testing  new,  miniature  technology  aboard 
its  three  microsatellites,  each  of  the  25  kg  class  [15].  The  satellites  were  delivered  into  a 
low-Earth  polar  orbit  with  an  altitude  of  4500  km.  aboard  the  same  launch  vehicle,  made 
possible  by  their  small  size.  The  satellites  were  then  configured  into  an  in-plane  forma¬ 
tion  as  they  made  measurements  of  the  Earth’s  magnetosphere,  proving  that  the  small, 
inexpensive  spacecraft  were  highly  capable  of  performing  the  mission. 
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The  Air  Force  continues  to  push  forward  with  formation  flying  experiments  as  well.  In 
1995,  the  Air  Force  New  World  Vistas  Space  Technology  panel  saw  the  potential  presented 
by  formation  flying,  advocating  the  exploration  of  “the  technical  challenges  and  benefits 
of  replacing  large  single  satellites  with  formations  of  microsatellites  to  perform  the  same 
mission”  [12].  This  initiative  led  directly  to  the  Air  Force  Research  Laboratory  (AFRL) 
TechSat  21  flight  experiment,  designed  to  test  and  validate  formation  flying  concepts  using 
three  microsatellites  that  function  as  one  “virtual  satellite.”  Although  the  mission  was 
cancelled  prior  to  launch,  its  formulation  spurred  considerable  interest  in  formation  flight 
and  led  to  increased  research  in  autonomous  formation  control.  AFRL  continues  to  push 
forward  the  concepts  of  satellite  formations  through  in-house  research  and  by  sponsoring 
microsatellite  initiatives  [1]. 

1-4  Problem  Description 

Research  involving  satellite  formation  flying  continues  to  progress.  With  technology 
advancing  to  facilitate  the  design  and  implementation  of  smaller  satellites,  fuel  concerns 
begin  to  dominate  the  cost  and  feasibility  of  missions.  Current  control  costs  are  not 
realistic  to  maintain  a  long-term  formation.  Additionally,  the  need  for  extensive  fuel 
capacity  dominates  the  design  of  the  satellite,  compromising  the  conception  of  smaller, 
cheaper  systems.  As  a  result,  methods  for  reducing  the  fuel  budget  without  compromising 
the  life  of  the  mission  lie  at  the  center  of  discussion.  Atmospheric  drag  is  potentially  an 
effective  actuator  for  formation  control  and  carries  the  added  benefit  of  being  in  unlimited 
supply  to  all  near-Earth  orbiting  satellites.  Without  the  need  for  expendables,  the  cost 
of  developing  and  launching  the  satellite  is  greatly  reduced.  However,  like  most  control 
techniques,  drag  has  its  limits.  Subjecting  a  spacecraft  to  excess  drag  increases  the  rate 
of  decay  the  satellite  would  otherwise  experience.  Therefore,  it  is  important  to  limit  the 
amount  of  drag  acting  on  a  spacecraft,  especially  when  drag  is  used  as  an  actuator. 

This  thesis  will  apply  a  linear-quadratic  terminal  feedback  control  law  to  reconfigure 
a  two-satellite  formation  using  atmospheric  drag  actuated  control  while  minimizing  the 
change  in  energy  of  the  satellites.  Since  the  orbital  energy  of  a  spacecraft  decreases  with 
orbital  radius,  minimizing  the  energy  change  and,  as  a  direct  result,  the  amount  of  control 
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used  will  maximize  the  lifetime  of  the  formation.  The  formation  will  be  reconfigured  over 
a  user-defined  fixed  time  about  an  empty  reference  orbit.  Each  satellite  is  equipped  with 
drag  panels  that  can  be  rotated  from  0  degrees,  which  is  parallel  to  the  flow  and  results  in 
the  minimum  drag  acceleration,  to  90  degrees,  which  is  normal  to  the  flow  and  results  in  the 
maximum  drag  acceleration.  The  rotation  of  the  drag  panels  effectively  varies  the  presented 
area  of  the  satellite,  which  is  directly  proportional  to  the  ballistic  coefficient.  A  maximum 
differential  ballistic  coefficient  of  0.1  m2 /kg  is  used  for  the  study,  a  conservative  value  for 
the  control  input.  Reconfigurations  to  both  an  in-plane  and  an  elliptical  formation  will 
be  performed  from  three  initial  formations:  an  in-plane  formation,  an  elliptical  formation, 
and  a  post-payload  separation  configuration.  The  performance  of  each  reconfiguration  will 
be  evaluated  in  terms  of  altitude  lost  and  total  Av  required,  the  typical  fuel  budget  for 
orbital  transfers.  Moreover,  the  final  conditions  of  the  controlled  state  will  be  propagated 
forward  in  time  to  verify  the  terminal  constraints  of  the  control  law  were  met. 

The  primary  goal  of  the  study  is  to  show  that  a  satellite  formation  can  be  effectively 
controlled  using  atmospheric  drag  while  minimizing  the  loss  of  altitude  in  the  process.  As 
such,  several  assumptions  were  made  to  help  focus  the  study  in  this  respect: 

•  The  reference  position  is  following  a  circular,  equatorial  orbit  with  the  given  orbital 
radius; 

•  The  atmosphere  is  assumed  to  decay  exponentially  with  increased  altitude  and  rotates 
with  the  Earth; 

•  The  atmospheric  density  at  the  reference  altitude  is  constant  throughout  the  orbit 
and  remains  constant  over  the  duration  of  the  reconfiguration; 

•  The  molecular  collisions  transfer  all  energy  in  the  in-track  direction.  Therefore,  no 
lift  is  generated,  just  a  retarding  force  against  the  satellite’s  motion; 

•  The  satellites  have  identical  ballistic  coefficients  when  their  drag  panels  are  at  the 
same  orientation,  and  the  drag  panel  orientation  does  not  affect  the  attitude  of  the 
satellite; 

•  The  full  state  of  each  satellite  can  be  measured  for  the  feedback  control  law; 
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•  The  relative  velocities  of  the  satellite  do  not  contribute  to  the  drag  acceleration;  only 
the  orbital  velocity  of  the  circular  reference  orbit  is  considered  since  it  is  on  the  order 
of  104  larger  than  the  relative  velocities. 

1.5  Thesis  Outline 

The  dynamics  of  the  study  are  discussed  in  Chapter  II.  In  Section  2.1,  the  relative 
equations  of  motion  are  derived,  defining  the  position  and  motion  of  a  satellite  with  respect 
to  a  reference  position  following  a  circular  orbit.  Then,  the  solution  to  these  equations  of 
motion  are  developed  for  the  unforced  and  forced  motion  to  give  insight  into  the  natural 
dynamics  of  the  problem.  The  two  common  formation  design  schemes  are  defined,  along 
with  the  advantages  and  disadvantages  of  each  design.  The  effects  of  atmospheric  drag 
on  satellite  motion  is  outlined  in  Section  2.2,  with  a  brief  discussion  of  the  exponential 
atmospheric  density  model  used  for  the  study  and  the  derivation  of  the  drag  acceleration 
imparted  on  a  satellite.  The  orbital  energy  and  the  change  in  energy  due  to  drag  is  defined 
in  Section  2.3.  And  finally  in  Section  2.4,  the  falling  and  fixed  reference  frames  used  for 
this  study  are  presented. 

Chapter  III  outlines  the  control  law  used  for  the  reconfigurations.  In  Section  3.1, 
the  controllability  of  a  single  satellite  using  atmospheric  drag  actuated  control  is  used  to 
define  the  relevant  state  equations  of  a  single  satellite  for  this  study.  Since  the  cross-track 
direction  is  uncontrollable  using  drag,  only  the  radial  and  in-track  components  of  each 
satellite  are  needed  for  the  controller.  Also,  the  control  law  is  dependent  on  the  states 
of  both  satellites,  so  the  state  equations  of  each  satellite  are  combined  in  Section  3.2  to 
give  the  equations  used  by  the  controller.  The  optimal  control  law  is  derived  in  Section 
3.3  from  the  Euler-Lagrange  necessary  conditions,  with  the  performance  index  chosen  to 
penalize  the  loss  of  orbital  energy  and  to  weight  the  difference  between  the  desired  and 
actual  final  state.  Finally,  the  constraints  needed  for  the  control  law  to  effectively  drive 
the  satellites  into  a  final  in-plane  and  elliptical  formation  with  zero  differential  drift  are 
derived  in  Section  3.4. 

Chapter  IV  contains  the  results  of  the  simulation.  The  numerically  integrated  equa¬ 
tions  of  motion  are  first  verified  against  the  closed  form  solution  for  the  unforced  motion 
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and  a  linear  approximation  for  the  forced  motion  in  Section  4. 1 .  The  three  reconfigurations 
to  a  final  in-plane  formation  are  presented  in  Section  4.2,  and  reconfigurations  to  a  final  el¬ 
liptical  formation  are  presented  in  Section  4.3.  The  time  history  and  radial  versus  in-track 
position  plots  for  both  the  reconfiguration  and  post-reconfiguration  analysis  are  presented 
within  the  two  sections  to  give  insight  into  the  motion  of  the  satellites  and  to  ensure  the 
final  constraints  were  met  within  sufficient  bounds.  A  Av  budget  is  also  provided  for  each 
reconfiguration  as  a  means  of  comparison  for  conventionally  controlled  formations. 

Finally,  Chapter  V  contains  the  conclusions  of  the  study  and  recommendations  for 
future  work. 
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II.  Background 


2.1  Relative  Satellite  Motion 

Within  this  section,  the  relative  equations  of  motion  applied  to  this  study  will  be 
developed,  as  well  as  the  unforced  and  forced  solutions  to  these  equations.  Additionally, 
common  formation  types  and  designs  will  be  discussed. 

2.1.1  Relative  Equations  of  Motion.  To  derive  the  equations  of  motion,  a  cylin¬ 
drical  coordinate  frame  will  be  used.  As  shown  in  Figure  2.1,  the  radial  direction,  er,  is 
aligned  with  the  position  vector  of  the  reference  position  and  is  measured  radially  outward. 
The  in-track  direction,  eg,  is  perpendicular  to  er  in  the  direction  of  the  velocity  vector  of 
the  reference  position.  For  circular  orbits,  eg  is  measured  along  the  velocity  vector  of  the 
reference  position.  And  finally,  the  cross-track  direction,  ez,  is  normal  to  the  orbital  plane, 
completing  the  right-handed  set. 

Using  the  coordinate  system  defined  in  Figure  2.1,  the  position  vector  of  the  reference 
position,  rref,  can  be  written  as 


r  ref  =  r0er  (2.1) 

where  ro  is  the  orbital  radius  of  the  reference  position,  which  follows  a  circular  orbit.  Define 
the  position  of  a  satellite  relative  to  the  reference  position,  psat,  such  that 


psat  =  5rer  +  r059ee  +  dzez  (2.2) 

where  dr,  r^dO,  and  dz  are  the  radial,  in-track,  and  cross-track  displacements  from  the 
reference  position,  respectively.  Given  the  definition  of  the  relative  position,  the  absolute 
position  of  the  satellite,  rsat,  can  be  written  as 


rsat.  —  rref  T  Psat. 

=  (ro  +  dr)er  +  rodOeg  +  dzez 


(2.3) 
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Figure  2.1  Reference  Frame 

This  position  vector  is  written  in  terms  of  the  relative  frame,  which  rotates  with  the  circular 
reference  orbit  at  the  constant  rate  of  n,  the  mean  motion  of  the  reference  orbit.  So  the 
angular  velocity  and  acceleration  of  the  reference  frame  are  defined  as 


uei  =  nez  (2.4) 

uei  =  0  (2.5) 

(2.6) 

where  u)ei  is  zero  since  the  angular  velocity  vector  is  constant,  and  fi  is  the  gravitational 
parameter  of  Earth.  The  absolute  acceleration  of  the  satellite  with  respect  to  the  relative 
frame  is  then  defined  as 
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ed2 

^sat  —  , ,  o  J'sat  “1“  2tt? 


ed 
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where  sat  and  rsat  are  the  time  derivatives  of  the  satellite  position  vector  in  the 
relative  frame  only,  such  that 


5rer  +  rodOeg  +  dzez  (2-8) 

5rer  +  r^ddeg  +  dzez  (2-9) 

Substituting  the  expressions  from  Equations  2.4,  2.8,  and  2.9  into  Equation  2.7  and  equat¬ 
ing  the  cross  products  yields  the  acceleration  vector  of  the  satellite  in  terms  of  the  relative 
frame  components: 


7, r sat  — 

at 


y2 

dt 2 


I*  sat  — 


rsat  =  (Sr  —  2nroS8  —  n2?’o  —  n2dr)er  +  (r^SO  +  2  ndr  —  n2roS0)eg  +  dzez  (2.10) 


Assuming  two-body  motion,  another  expression  for  the  acceleration  vector  of  the 
satellite  can  be  written  as 


rSat  =  ~-^-rsat  +  f  (2.11) 

^  sat 

The  first  term  in  Equation  2.11  is  the  two-body  gravitational  acceleration  acting  on  the 
satellite.  The  vector  f,  which  can  be  written  in  component  form  such  that  f  =  frer  + 
fgeg  +  fzezi  contains  any  forces,  per  unit  mass,  other  than  the  two-body  gravitational 
force  acting  on  the  satellite.  Substituting  the  position  vector  of  the  satellite  from  Equation 
2.10  into  Equation  2.11  yields 


H[(r0  +  dr)er  +  r0ddeg  +  dzez] 
(tq  +  2r$dr  +  dr2  +  r^dO2  +  Sz2)3/2 


(2.12) 
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Assuming  small  relative  displacements  with  respect  to  the  orbital  radius,  ro,  the  denomi¬ 
nator  of  Equation  2.12  can  be  expanded  using  the  binomial  theorem,  resulting  in 

(vq  +  2r0Sr  +  dr2  +  r%592  +  8z2)  '^2  =  r^3  1  ~  +  •  •  •  ^  +  0(82)  (2.13) 

Substituting  this  result  into  Equation  2.12  and  neglecting  terms  of  order  0(8 2)  and  higher 
yields 

r sat  =  --^[(ro  -  2 5r)er  + r059ee  +  8zez }  +  f  (2.14) 

ro 

Combining  Equations  2.10  and  2.14  and  recalling  the  definition  of  the  mean  motion  results 
in  the  following  three  equations: 


5r  = 

2nro89  +  3  n2dr  +  fr 

(2.15) 

ro59  = 

-2  n8r  +  fe 

(2.16) 

8z  = 

- r?8z  +  fz 

(2.17) 

These  three  equations  are  commonly  known  as  Hill’s  equations  or  the  Clohessy- Wiltshire 
(CW)  equations  and  are  used  to  describe  the  relative  motion  of  two  satellites  when  the 
reference  satellite  is  following  a  circular  orbit  [25].  The  double  name  credits  both  parties 
involved  in  the  development  of  the  equations.  Hill  first  studied  the  equations  in  1878,  but 
Clohessy  and  Wiltshire  redeveloped  the  equations  in  1960  for  satellite  rendezvous  [24,  4]. 

In  forming  these  equations,  several  assumptions  were  made.  One  assumption  is  that 
the  orbital  mean  motion  n  is  constant.  This  assumption  is  valid  for  circular  orbits,  but  can 
be  extended  to  nearly  circular  orbits  within  bounds.  Another  assumption  was  made  when 
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Figure  2.2  Illustration  of  in-track  measurement 


substituting  the  result  of  the  binomial  expansion  from  Equation  2.13  into  the  two-body 
gravitational  acceleration  in  Equation  2.14.  Here,  terms  of  order  0(S2)  and  higher  were 
neglected,  an  assumption  valid  for  small  coordinates  only.  However,  the  final  equations 
are  not  limited  in  the  tq59  direction  since  this  quantity  does  not  appear  explicitly  in 
the  equations  above  and  is  measured  circumferentially,  as  illustrated  in  Figure  2.2  [7]. 
This  extended  range  of  validity  in  the  in-track  direction  is  a  benefit  of  using  cylindrical 
coordinates. 

2.1.2  Force  Free  Solution.  To  further  understand  the  relative  motion  of  the 
satellite,  the  unforced  solution  to  the  CW  equations  will  be  developed  following  the  same 
process  used  by  Wiesel  [25].  For  this  case,  f  =  0,  resulting  in  the  following  equations: 


Sr  = 

2nroS9  +  3  n2Sr 

(2.18) 

r(j80  = 

—2  nSr 

(2.19) 

Sz  = 

—n2Sz 

(2.20) 
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As  can  be  seen  from  the  equations  above,  the  8z  equation  is  decoupled  from  the  in-plane 
equations  of  motion  and  represents  the  motion  of  a  harmonic  oscillator.  So,  given  the 
initial  conditions  5z(t  =  0)  =  Szo  and  5z(t  =  0)  =  Szq,  the  solution  is  simply 

&  z 

Sz(t)  =  Szo  cos  nt  H - sin nt  (2.21) 

n 

Integrating  Equation  2.19  with  initial  conditions  59(t  =  0)  =  56 o  and  Sr(t  =  0)  =  Sro  and 
solving  for  56  yields 


S6(t)  =  560  +  —(5r0-5r)  (2.22) 

Substituting  this  equation  into  the  5r  equation  results  in 

Sr  +  n25r  =  2nro59o  +  4n2(5ro  (2.23) 

This  equation  is  a  simple  harmonic  oscillator  with  a  forcing  function.  Its  solution  contains 
both  a  homogeneous  solution  and  a  particular  solution  such  that 

5r(t)  =  5rh  +  5rp  (2.24) 

where  the  particular  solution  is  just 


5rp  =  -r0590  +  4 5r0 
n 

Solving  for  the  homogeneous  solution  with  initial  conditions 


(2.25) 


Oq 

c-t- 

II 

=  5r0  -  5rp  =  ~-r0560  -  35r0 
n 

(2.26) 

Srh(t  =  0) 

=  5r0 

(2.27) 
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results  in  the  complete  solution 


5r(t)  =  —  ( — ro66o  +  35?’o  )  cos  nt  H - —  sin  nt  H — ro50o  +  45ro  (2.28) 

\n  J  n  n 

Substituting  this  solution  into  Equation  2.22  gives 

•  •  (St?  {  •  6r?  \ 

56(t)  =  —3<50o - 5ro  +  (  4(50o  H - <5ro  )  cos  nt - sin  nt  (2.29) 

ro  V  r0  J  r0 

Integrating  both  sides  of  Equation  2.29,  applying  initial  conditions,  and  solving  for  5Q(t ) 
yields 

(50(f)  =  56q  —  [350o  4 - (5r0  |  t  +  ( — 56q  H - 5?’o  )  sin  rat  +  0  cosnf - 5ro  (2.30) 

V  r0  )  \n  r0  )  nr0  nr0 

Equations  2.21,  2.28,  and  2.30  form  the  position  solution  in  terms  of  an  initial  state  and 
final  time.  To  find  the  velocity  solution,  simply  differentiate  the  position  solution,  yielding 

5r(t)  =  (2ro<50o  +  3n<5?’o)  sin  nt  +  <5?’o  cosnf  (2.31) 


5z{t )  =  —  fco^  sin  nt  +  Szq  cos  nt  (2.32) 

and  where  50(f)  is  given  in  Equation  2.29.  Equations  2.29,  2.31,  and  2.32  represent  the 
velocity  solution  to  the  unforced  equations  of  motion.  Combining  the  position  and  velocity 
solutions,  the  overall  solution  can  be  put  into  state  variable  form,  producing 
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4  —  3  cos  nt 

0 

0 

1  sin  nt 

n 

4(1  —  cos  nt) 
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4  (cos  nt  —  1) 
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0 

0 

cos  nt 
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4  sin  nt 

Sz0 

Sr 

3n  sin 
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0 

cos  nt 

2  sin  nt 
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Sf0 
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6n(cos  nt  —  1) 

0 
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—2  sin  nt 

—3  +  4  cos  nt 
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Sz 
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0 
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0 

cos  nt 

O 

(2.33) 


Several  trends  are  evident  from  this  solution.  First,  the  cross-track  components,  Sz 
and  Sz,  simply  oscillate  with  the  period  of  the  reference  orbit.  This  oscillation  results  in 
a  slight  inclination  for  the  satellite,  projecting  it  into  a  separate  orbital  plane  from  the 
reference  orbit.  Also  notice  that  the  cross-track  components  are  completely  decoupled 
from  the  in-plane  components,  making  all  motion  in  the  cross-track  direction  independent 
of  any  in-plane  motion.  The  radial  components,  Sr  and  Sr,  contain  the  same  oscillatory 
terms  as  the  cross-track  components  that  move  with  the  period  of  the  reference  orbit. 
This  motion  creates  a  slight  eccentricity  in  the  satellite’s  orbit  as  the  radius  increases  and 
decreases.  Finally,  the  in-track  components,  r^SO  and  r^SO,  contain  the  same  oscillations, 
but  ro59  also  includes  a  term  that  grows  with  time.  This  secular  term  causes  an  in-track 
drift  between  the  reference  position  and  the  satellite.  To  maintain  a  formation  for  an 
extended  period  of  time,  this  secular  motion  must  be  suppressed. 

2.1.3  Forced  Solution.  Unlike  the  force  free  solution  presented  in  Section  2.1.2, 
a  nice,  closed-form  solution  does  not  exist  for  the  forced  CW  equations.  However,  one  can 
be  constructed  that  uses  parameters  similar  to  the  force  free  solution.  Following  the  same 
steps  used  by  Meirovitch  to  solve  a  nonhomogeneous  linear  system  of  equations  [13],  define 
the  forced  equations  of  motion  from  Equations  2.15-2.17  in  state  vector  form  such  that 

^x(t)  =  Ax(t)  +  i(t)  (2.34) 

where  x,  A,  and  f  are  defined  as 
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Sr  tqSO  Sz  Sr  tqSO 
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f 
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0  0  0  0  0  1 

3 n2  00  0  2n  0 

0  0  0  —2 n  0  0 

0  0  —  n2  0  0  0 


0  0  0 


fr  h  fz 


T 


(2.35) 


(2.36) 


(2.37) 


Introduce  the  state  transition  matrix,  d>(t),  which  satisfies  the  equations 


$(to)  =  /  (2-38) 

Jtm  =  A$(t)  (2.39) 

i-^-1^)  =  -Q-^A  (2.40) 

where  /  is  the  identity  matrix.  Note  the  inverse  of  the  state  transition  matrix  is  needed  in 
Equation  2.40.  This  inverse  exists  for  all  time  by  definition  of  the  Jacobi-Liouville  formula 
and  Equation  2.38  [13].  Premultiplying  Equation  2.34  by  <b~1(f)  and  postmultiplying 
Equation  2.40  by  x(t),  the  two  equations  can  be  added  together  to  yield 


|  [*-'(t)x(()]  =4>-1(«)f(*) 


(2.41) 
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Integrating  Equation  2.41  with  respect  to  time  and  premultiplying  by  4>(t)  yields  the 
solution  to  the  forced  equations  of  motion: 


x(t) 


<3?(t)x(f0)  +  4>(f)  [  $  1(r)f(r)dr 
JtQ 


where 


(2.42) 


m 


4  —  3  cos  nt  0 
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(2.43) 


2.1.4  Formation  Design.  The  design  of  a  satellite  formation  is  critical  to  the 
control  costs  to  maintain  the  formation.  Two  common  formation  schemes  exist:  a  leader- 
follower  scheme  and  an  empty  reference  position  formation,  both  of  which  were  alluded  to 
in  Section  1.1. 

The  most  common  design  is  the  leader-follower  scheme.  This  two-satellite  formation 
consists  of  the  leader  satellite,  which  is  usually  at  the  center  of  the  reference  frame,  and 
a  follower  satellite,  which  is  positioned  relative  to  the  leader.  The  leader  usually  takes 
a  passive  role,  meaning  it  is  non-maneuvering.  Instead,  it  flies  in  a  circular  orbit,  and 
the  follower  satellite  is  controlled  around  it,  taking  on  the  active  role.  This  scheme  is 
perhaps  the  most  common  due  to  its  similarities  to  orbital  rendezvous,  an  older  concept  in 
astro  dynamics.  For  rendezvous,  the  station  is  usually  in  a  fixed  orbit  and  takes  no  action 
in  docking  with  the  approaching  satellite.  And  since  the  CW  equations  were  developed 
for  rendezvous,  the  dynamics  translate  well  for  the  leader-follower  scheme,  although  not 
completely.  Since  rendezvous  is  usually  performed  over  a  short  time  period,  perturbing 
effects  can  be  ignored  with  confidence  since  they  usually  require  an  extended  period  of  time 
to  noticeably  impact  satellite  motion.  Also,  rendezvous  usually  involves  a  considerable 
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fuel  budget,  allowing  the  cancellation  of  any  disturbances  through  increased  control.  But 
when  describing  a  satellite  formation  that  must  maintain  its  relative  positioning  for  several 
months  to  years  at  a  time  with  a  limited  fuel  budget,  neglecting  these  perturbations  could 
be  costly,  especially  with  the  design  of  the  leader-follower  scheme.  As  is  evident  from 
the  design,  the  two  satellites  are  flying  in  separate  orbits,  allowing  any  orbital  element 
dependent  perturbations  to  pull  the  formation  apart.  One  such  disturbance  that  affects 
low-Earth  orbit  satellites  is  the  J2  perturbation,  the  geopotential  zonal  harmonic  associated 
with  the  oblateness  of  the  Earth  [26].  The  J2  perturbation  places  a  torque  on  the  orbit 
that  causes  the  orbital  plane  to  precess  and  the  argument  of  perigee  to  rotate  at  rates  of 
[25] 


n 


2 a2  (1  -  e2)2 


cosi 


(2.44) 


Cj 


3nJ2-Re 
2a?  (1  —  e2)2 


(5  .  2. 

y  2  1 
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where  Q  is  the  right  ascension  of  the  ascending  node,  co  is  the  argument  of  perigee,  J2  is  a 
dimensionless  number  with  the  value  of  0.001082  that  characterizes  the  Earth’s  oblateness, 
Re  is  the  radius  of  the  Earth,  a  is  the  semimajor  axis,  e  is  the  eccentricity,  and  i  is  the 
orbital  inclination.  The  J2  perturbation  tends  to  pull  satellite  formations  apart  if  the  two 
satellites  have  different  inclinations  or  eccentricities,  which  is  the  case  for  most  leader- 
follower  schemes.  As  an  example,  consider  an  elliptical  formation  with  the  leader- follower 
arrangement  such  that  the  leader  is  following  a  circular  orbit  and  the  follower  is  in  a 
slightly  eccentric  orbit,  allowing  it  to  orbit  the  leader  in  the  relative  frame.  Because  of 
their  eccentricity  difference,  the  two  orbital  planes  of  the  satellites  will  precess  at  slightly 
different  rates,  causing  the  satellites  to  drift  apart.  This  drift  rate  is  increased  if  a  cross¬ 
track  motion  is  also  applied  to  the  follower  since  the  inclinations  will  differ.  Therefore,  the 
leader-follower  formation  is  intrinsically  unstable,  requiring  formation  keeping  maneuvers 
to  maintain  the  formation. 
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The  empty  reference  position  design  is  another  formation  type  that  is  sometimes 
employed.  For  this  design,  the  reference  position  is  left  empty,  meaning  there  is  not  a 
satellite  at  the  center  of  the  formation.  Instead,  the  satellites  are  positioned  around  a 
defined  reference  point  centered  on  the  formation.  This  formation  design  is  more  stable 
than  the  leader-follower  scheme  since  the  satellites  usually  have  the  same  semimajor  axis, 
inclination,  and  eccentricity.  For  example,  consider  the  same  elliptical  formation  described 
for  the  leader-follower  scheme,  but  using  the  empty  reference  position  instead.  To  place 
the  two  satellites  in  the  correct  formation,  their  semimajor  axes,  eccentricities,  and  incli¬ 
nations  must  be  equal,  with  their  arguments  of  perigee  determining  their  placement  in  the 
formation.  As  a  result,  their  orbital  planes  and  arguments  of  perigee  would  rotate  at  the 
same  rate,  keeping  the  formation  intact.  The  reference  position  would  merely  rotate  at 
this  rate  as  well  to  maintain  the  average  position  of  the  formation.  As  long  as  the  overall 
drift  of  the  formation  is  accounted  for  in  the  mission  design,  the  empty  reference  position 
formation  scheme  provides  a  stable  formation  that  does  not  require  considerable  formation 
keeping  costs. 

2.2  Atmospheric  Drag 

The  effects  of  the  J2  zonal  harmonic  on  satellite  motion  was  discussed  in  the  previous 
section,  but  the  gravitational  forces  are  not  the  only  forces  affecting  satellite  motion.  The 
atmosphere,  although  extremely  thin  at  orbital  altitudes,  can  have  a  substantial  impact 
on  low-Earth  orbit  (LEO)  satellites.  In  fact,  besides  gravitational  forces,  atmospheric  drag 
is  the  most  influential  force  acting  on  LEO  satellites,  retarding  their  motion  and  leading 
to  eventual  decay  into  the  denser  layers  of  the  atmosphere  [9].  Within  this  section,  the 
atmospheric  model  used  for  this  study  will  be  introduced,  and  the  effect  of  drag  on  satellite 
motion  will  be  derived.  Finally,  the  impact  of  differential  drag  on  satellites  in  formation 
will  be  discussed. 

2.2.1  Atmospheric  Model.  Perhaps  the  most  difficult  aspect  of  the  near-Earth 
space  environment  to  model  is  the  atmosphere.  The  atmosphere  is  a  highly  dynamic 
medium  that  can  vary  significantly  at  a  given  altitude  over  time.  Certain  variations  are 
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predictable,  such  as  the  equatorial  bulge  and  the  diurnal  bulge  [26].  However,  atmospheric 
variations  due  to  solar  radiation  and  the  Earth’s  magnetic  field  can  be  highly  unpredictable, 
making  atmospheric  conditions  difficult  to  predict  into  the  future  [17,  26].  Therefore,  most 
models  are  simply  approximations  of  the  atmosphere  under  normal  conditions. 

Since  the  purpose  of  this  study  is  to  prove  the  feasibility  of  controlling  a  satellite 
formation  through  the  use  of  drag,  a  comprehensive  atmospheric  model  complete  with  time 
dependent  variations  is  not  needed.  Instead,  a  simple  model  for  estimating  the  atmospheric 
density  at  a  given  altitude  is  sufficient.  Therefore,  an  exponential  atmospheric  model  is 
used  for  this  study.  This  model  assumes  a  static  atmosphere  with  a  spherically  symmetrical 
distribution  of  particles  [22],  The  atmospheric  density,  p,  decreases  exponentially  with 
increasing  altitude,  varying  according  to  the  following  relationship: 


P  =  Poexp 


hellp  ho 

H 


(2.46) 


Here,  po  is  the  reference  density  at  the  reference  altitude  ho,  heup  is  the  altitude  measured 
from  the  Earth’s  surface,  and  H  is  the  scale  height  at  the  reference  altitude  [22].  Values 
for  these  definitions  are  given  in  Table  2.1. 


2.2.2  Acceleration  Due  to  Drag.  The  atmosphere  at  orbital  altitudes  consists  of 
diffuse  molecules  that  have  very  little,  if  any,  interaction  with  one  another.  This  region, 
known  as  the  free  molecular  flow  regime,  consists  of  molecules  that  no  longer  behave  like 
a  fluid,  but  instead  behave  as  individual  particles  [17].  At  orbital  speeds,  the  individual 
molecular  impacts  on  the  spacecraft,  when  taken  over  an  extended  period  of  time,  can 
greatly  affect  the  motion  of  the  satellite.  When  the  molecules  impact  and  bounce  off  the 
spacecraft,  linear  momentum  is  transferred  through  the  inelastic  collisions  proportional 
to  the  velocity  difference  between  the  molecules  and  the  spacecraft  [26].  Since  force  is 
the  rate  of  change  of  linear  momentum,  these  individual  impacts  transfer  a  drag  force 
to  the  spacecraft  that  imparts  an  acceleration  in  the  direction  opposite  the  spacecraft’s 
motion  proportional  to  Apv2 / 2.  Here,  A  is  the  presented  area  of  the  satellite,  p  is  the 
atmospheric  density  at  the  given  orbital  altitude,  and  v  is  the  satellite  velocity  with  respect 
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Table  2.1  Exponential  Atmospheric  Model  [22] 


Altitude 

hellp 

(km) 

Base 

Altitude 

ho 

Nominal 

Density 

p0(kg/m3) 

Scale 

Height 

H(km) 

0-25 

0 

1.225 

7.249 

25-30 

25 

3.899  x  10~2 

6.349 

30-40 

30 

1.774  x  10~2 

6.682 

40-50 

40 

3.972  x  10~3 

7.554 

50-60 

50 

1.057  x  10~3 

8.382 

60-70 

60 

3.206  x  10-4 

7.714 

70-80 

70 

8.770  x  10“5 

6.549 

80-90 

80 

1.905  x  10“5 

5.799 

90-100 

90 

3.396  x  10~6 

5.382 

100-110 

100 

5.297  x  10~7 

5.877 

110-120 

110 

9.661  x  10”8 

7.263 

120-130 

120 

2.438  x  10~8 

9.473 

130-140 

130 

8.484  x  10~9 

12.636 

140-150 

140 

3.845  x  10-9 

16.149 

150-180 

150 

2.070  x  10“9 

22.523 

180-200 

180 

5.464  x  10~10 

29.740 

200-250 

200 

2.789  x  10~10 

37.105 

250-300 

250 

7.248  x  10~n 

45.546 

300-350 

300 

2.418  x  10~n 

53.628 

350-400 

350 

9.518  x  10~12 

53.298 

400-450 

400 

3.725  x  10“12 

58.515 

450-500 

450 

1.585  x  10“12 

60.828 

500-600 

500 

6.967  x  10“13 

63.822 

600-700 

600 

1.454  x  10”13 

71.835 

700-800 

700 

3.614  x  10"14 

88.667 

800-900 

800 

1.170  x  10"14 

124.64 

900-1000 

900 

5.245  x  10-15 

181.05 

1000- 

1000 

3.019  x  10”15 

268.00 
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to  the  atmosphere.  The  shape  of  the  satellite  also  affects  the  transfer  of  momentum  to 
the  spacecraft.  As  a  result,  the  resultant  drag  force  is  also  directly  proportional  to  the 
coefficient  of  drag,  Cd-  This  aerodynamic  coefficient  is  a  function  of  body  shape  and 
the  Mach  number,  along  with  other  various  parameters,  and  also  accounts  for  oblique 
molecular  reflections  [17,  26].  Thus,  the  acceleration  due  to  drag,  a^,  is  written  as 


a  d  =  ~ 


1  CdA 

2  m 


-pvw 


(2.47) 


where  m  is  the  satellite  mass.  Assuming  a  circular,  equatorial  orbit  with  an  atmosphere 
that  rotates  with  the  Earth,  the  satellite  velocity  vector  with  respect  to  the  atmosphere, 
v,  is  defined  as  [26] 


v  =  vin  -  uje  x  r  (2.48) 

where  vm  is  the  inertial  velocity  of  the  satellite,  uie  is  the  rotational  velocity  of  the  Earth, 
and  r  is  the  inertial  satellite  position  vector. 

The  drag  coefficient,  presented  area,  and  mass  may  not  be  separately  determinable, 
so  these  three  quantities  are  usually  grouped  into  a  single  quantity  called  the  ballistic 
coefficient,  B*,  which  is  defined  as  [26] 


B,  =  CnA 

m 


(2.49) 


From  this  definition,  it  can  be  seen  that  increasing  the  ballistic  coefficient  increases  the 
amount  of  drag  acting  on  the  satellite.  Since  the  drag  coefficient  is  relatively  fixed,  the 
ballistic  coefficient  can  change  only  if  the  presented  area  of  the  satellite  or  the  satellite  mass 
is  changed.  Changing  the  satellite  mass  is  an  irreversible  process  since  it  usually  involves 
the  bleeding  of  fuel,  so  a  change  in  mass  is  not  an  ideal  method  for  altering  the  ballistic 
coefficient  for  control  purposes.  However,  the  presented  area  of  the  satellite  can  easily 
change  by  a  simple  rotation  of  the  satellite  itself  or  through  the  rotation  of  either  drag 
panels  or  solar  panels  attached  to  the  satellite,  making  it  an  ideal  technique  for  control. 
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Note  that  several  conventions  exist  for  defining  the  ballistic  coefficient,  one  being 
the  inverse  of  the  definition  above  such  that  units  of  pressure  result  [22],  However,  the 
definition  in  Equation  2.49  is  an  interpretation  developed  in  Wiesel’s  texts  and  is  the 
accepted  definition  for  this  study  [25,  26]. 

2.2.3  Differential  Drag.  When  controlling  a  satellite  formation  with  drag,  one 
satellite  is  purposely  subjected  to  increased  drag  to  change  its  position  with  respect  to 
the  other.  The  drag  force  difference  between  the  two  satellites  is  called  differential  drag, 
and  the  force  difference  results  in  a  differential  acceleration  that  can  be  either  negative  or 
positive  in  the  relative  frame,  although  the  drag  acting  on  the  satellite  is  always  directed 
against  the  satellite’s  velocity  vector. 

To  help  understand  the  concept  of  differential  drag,  consider  two  identical  satellites, 
both  equipped  with  drag  panels,  flying  in  the  same  orbit  but  separated  by  true  anomaly 
such  that  satellite  1  is  slightly  ahead  of  satellite  2  in  an  in-plane  formation.  Assume  satellite 
1  is  at  the  center  of  the  relative  frame.  If  both  satellites  have  their  drag  panels  oriented  at 
zero  degrees  such  that  the  drag  force  acting  on  each  satellite  due  to  the  panels  is  zero,  the 
satellites  experience  zero  differential  drag.  They  will  neither  move  toward  nor  away  from 
one  another.  If  the  panels  on  satellite  1  were  then  rotated  to  90  degrees  such  that  they  were 
normal  to  the  velocity  vector,  satellite  1  would  experience  a  negative  acceleration  in  the 
absolute  frame,  causing  satellite  2  to  move  toward  it.  As  a  result,  satellite  2  experiences 
a  positive  acceleration  in  the  relative  frame.  If  the  panels  on  satellite  2  were  then  rotated 
such  that  both  satellites  had  their  panels  oriented  at  90  degrees,  the  differential  drag  would 
then  be  zero  since  both  satellites  are  experiencing  the  same,  albeit  larger,  drag  force.  If  the 
panels  on  satellite  1  were  then  rotated  back  to  zero  degrees  while  the  panels  on  satellite  2 
stayed  at  90  degrees,  satellite  2  would  then  experience  a  negative  differential  acceleration 
since  it  would  then  move  away  from  satellite  1. 

Now  consider  satellite  1  with  its  panels  rotated  at  90  degrees  and  satellite  2  with 
its  panels  oriented  at  45  degrees.  The  differential  acceleration  experienced  with  this  con¬ 
figuration  is  the  same  differential  acceleration  experienced  when  the  panels  on  satellite  1 
are  oriented  at  45  degrees  and  the  panels  on  satellite  2  are  at  zero  degrees,  although  the 
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first  configuration  is  experiencing  double  the  drag.  This  increased  drag  does  nothing  but 
increase  the  rate  of  decay  of  the  formation  altitude,  decreasing  the  overall  energy  of  the 
formation.  As  a  result,  to  minimize  the  loss  of  altitude  and  the  loss  of  orbital  energy,  only 
one  satellite  should  have  its  panels  deployed  at  any  given  time.  Any  control  law  using 
differential  drag  as  an  actuator  should  account  for  this  fact. 

2. 3  Energy  Function 

When  controlling  a  spacecraft  with  atmospheric  drag,  the  satellite  is  purposely  sub¬ 
jected  to  increased  drag  to  achieve  the  desired  result.  Since  drag  acts  against  the  motion 
of  the  vehicle,  it  instantaneously  decreases  the  velocity  of  the  satellite,  which  causes  the 
orbital  radius  to  decrease.  Although  the  instantaneous  velocity  of  the  satellite  decreases, 
the  average  velocity  of  the  satellite  increases  as  a  result  of  the  reduction  in  orbital  size.  The 
continuous  decrease  in  orbital  size  will  shorten  the  lifetime  of  the  satellite  since,  without 
any  means  to  raise  the  orbit,  the  satellite  will  continue  to  decay  until  reentering  the  denser 
layers  of  the  atmosphere,  falling  out  of  orbit.  Therefore,  it  is  desirable  to  minimize  this 
rate  of  degradation  while  achieving  the  desired  result  of  the  control.  One  way  to  minimize 
the  degradation  is  to  minimize  the  loss  of  orbital  energy  of  the  satellite.  The  orbital  energy, 
£,  is  defined  as 


v2  fj,  n 
2  r  2  a 


(2.50) 


Here,  v  and  r  are  the  orbital  velocity  and  radius,  respectively,  and  a  is  the  semi-major 
axis.  As  is  evident  from  Equation  2.50,  reducing  the  orbital  size  of  the  satellite  reduces 
the  overall  energy.  And  since  drag  acts  to  shrink  the  orbital  size,  minimizing  the  change 
in  orbital  energy  for  this  case  is  the  same  as  minimizing  the  loss  of  orbital  energy.  The 
rate  of  change  of  orbital  energy  due  to  atmospheric  drag  can  be  found  by  differentiating 
Equation  2.50: 


d£  n 

—  =  vv  +  -~r 
at  rz 


(2.51) 
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Since  drag  acts  against  the  velocity  vector,  the  instantaneous  change  in  energy  due  to 
atmospheric  drag  becomes 


^  =  ~vad  (2.52) 

where  v  is  the  satellite  velocity  with  respect  to  the  atmosphere,  and  ad  is  the  acceleration 
due  to  drag. 

2-4  The  Relative  Frame 

The  relative  frame  used  for  this  study  will  be  slightly  different  from  the  reference 
frame  discussed  in  Section  2.1.1  and  shown  in  Figure  2.1.  The  reference  frame  discussed 
earlier  assumes  the  reference  position  is  in  a  fixed,  circular  orbit  such  that  its  orbital  radius 
and  velocity  are  constant.  However,  for  this  study,  atmospheric  drag  is  used  as  the  only 
means  for  control.  As  a  result,  the  orbital  radius  of  the  reference  position  will  decrease 
with  increased  control  input,  resulting  in  a  falling  reference  frame.  The  reference  position 
is  defined  as  the  average  position  of  the  two  satellites  in  the  formation,  so  it  will  fall  at  the 
instantaneous  average  rate  of  descent  of  the  satellites.  For  example,  suppose  one  satellite 
has  its  panels  oriented  such  that  the  drag  force  from  its  panels  is  zero  while  the  other 
satellite  has  its  panels  deployed  to  result  in  a  drag  acceleration,  causing  it  to  fall  at  a  rate 
of  r.  The  rate  of  change  of  the  reference  position  is  defined  as  r/2  such  that  the  reference 
position  remains  centered  on  the  formation. 

As  a  result  of  this  falling  frame  and  as  discussed  in  Section  2.2.3,  some  control  inputs 
appear  to  result  in  a  positive  acceleration  when  plotted  in  the  frame.  This  perception  is 
from  the  frame’s  motion.  A  reconfiguration  with  an  initial  20  m  displacement  in  the  radial 
direction  to  a  final  in-plane  formation  with  a  separation  distance  of  100  m.  is  shown  in  two 
different  relative  frames:  a  frame  fixed  upon  the  original  reference  position  in  Figure  2.3 
and  the  falling  frame  in  Figure  2.4.  The  fixed  frame  clearly  shows  that  the  final  reference 
position  fell  62  m  in  the  radial  direction  and  advanced  3400  m  in  the  in-track  direction  as  a 
result  of  the  reconfiguration.  But  the  position  of  the  satellites  with  respect  to  the  reference 
position  is  difficult  to  determine.  The  falling  frame  shows  the  more  common  relative  frame 


2-18 


In-track  (m) 

Figure  2.3  Plot  of  Reconfiguration  in  the  Fixed  Frame 

plot  with  a  reference  position  centered  on  the  formation  and  fixed  in  the  frame.  This 
plot  illustrates  the  appearance  of  a  positive  acceleration.  When  the  leading  satellite  has  a 
larger  drag  acceleration  than  the  trailing  satellite,  the  trailing  satellite  in  turn  experiences 
a  positive  acceleration  in  the  relative  frame.  In  reality,  the  acceleration  is  still  negative, 
but  the  drag  on  the  trailing  satellite  is  less  than  the  drag  on  the  leading  satellite. 

Also,  by  design  of  the  relative  frame,  the  forces  acting  on  the  two  satellites  in  the 
frame  are  differential  forces  only.  Equal  forces  acting  on  the  two  satellites  will  result  in  the 
same  motion  for  each  such  that  no  motion  results  when  viewed  in  the  relative  frame.  As 
an  example,  consider  two  satellites  flying  in  an  in-plane  formation  such  that  they  occupy 
the  same  circular  orbit  in  the  absolute  frame  but  are  separated  by  true  anomaly.  The  two 
satellites  will  continue  in  this  formation  as  long  as  all  forces  acting  on  the  two  satellites 
are  equivalent.  If  one  satellite  has  a  larger  drag  acceleration  than  the  other,  the  two  will 
separate  indefinitely  unless  a  control  input  is  applied. 

To  show  that  equivalent  forces  cancel  in  the  relative  frame,  consider  two  identical 
satellites  that  each  have  the  same  nominal  ballistic  coefficient  Bq  when  their  solar  panels 
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In-track  (m) 


Figure  2.4  Plot  of  Reconfiguration  in  the  Falling  Frame 

are  at  a  zero  degree  orientation  such  that  the  drag  force  due  to  the  panels  is  negligible. 
The  equations  of  motion  of  the  two  satellites  in  the  absolute  frame  then  become 

ri  =  --^ri  -  ^  (Bq  +  AB^pvv  (2.53) 

ri  l 

r2  =  -4r2  -  I,  (Bo  +  AB2)  pw  (2.54) 

r2  Z 

where  A B\  and  A are  the  ballistic  coefficients  due  to  non-zero  solar  panel  orientations 
for  satellite  1  and  satellite  2,  respectively.  Substituting  Equation  2.3  for  the  position 
vectors  of  each  satellite  and  expanding  the  denominator  using  the  binomial  theorem,  the 
resulting  in-track  accelerations  can  be  subtracted  in  the  relative  frame,  yielding 

rQSh  ~  ro^2  =  —  2n  (Sri  —  <5r2)  —  -  (A B\  —  A R|)  pv2  (2.55) 
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As  can  be  seen  from  Equation  2.55,  the  nominal  ballistic  coefficients  of  the  two  satellites 
cancel  one  another  in  the  relative  frame  when  they  are  equal,  and  only  the  differential 
ballistic  coefficient  due  to  the  orientation  of  the  panels  affects  the  motion  of  the  formation 
such  that  the  only  external  acceleration  seen  in  the  relative  frame  is  the  drag  acceleration 

ad  =  ^(AB*1-AB$)pv2  (2.56) 

This  same  argument  can  be  extended  to  the  gravitational  perturbations  acting  on 
the  formation.  As  long  as  the  two  satellites  are  identical  and  are  in  a  tight  formation,  the 
differential  forces  acting  upon  the  formation  will  be  negligible,  allowing  the  two  satellites 
to  remain  in  formation  with  little  disturbance. 


2-21 


III.  Control  Law 


When  considering  the  optimal  control  of  a  linear  system  using  quadratic  performance  crite¬ 
ria,  two  linear  quadratic  control  laws  come  to  mind:  a  regulator  and  a  terminal  controller. 
A  regulator  is  designed  to  maintain  a  system  within  an  acceptable  deviation  of  some  refer¬ 
ence  condition  using  acceptable  amounts  of  control,  while  a  terminal  controller  is  designed 
to  drive  a  system  from  some  initial  state  to  the  desired  final  state  over  a  finite  time  period, 
exhibiting  acceptable  behavior  along  the  way  [3].  The  purpose  of  this  study  is  to  opti¬ 
mally  reconfigure  a  two-satellite  formation,  driving  the  satellites  from  some  initial  state 
to  a  final  state  dependent  upon  the  desired  formation  while  minimizing  the  loss  of  orbital 
energy  of  the  formation.  Since  the  goal  is  to  reconfigure  the  formation,  not  to  perform 
station  keeping  of  a  current  formation,  a  terminal  controller  was  selected.  This  is  not 
to  say  that  a  linear  quadratic  regulator  (LQR)  control  law  cannot  be  used  to  reconfigure 
satellite  formations.  In  fact,  as  mentioned  in  Chapter  I,  several  instances  exist  where  a 
LQR  control  law  was  used  for  reconfigurations.  These  control  laws  usually  weight  both 
the  state  error  and  the  control  usage  within  the  performance  index  at  every  point  in  time, 
increasing  the  weight  on  the  state  error  matrix  to  decrease  the  time  required  to  perform 
the  reconfiguration,  or  increasing  the  weight  of  the  control  usage  penalty  matrix  to  reduce 
the  amount  of  control  used  to  drive  the  state  error  to  zero.  This  study  is  not  concerned 
with  a  time  optimal  reconfiguration.  Instead,  the  goal  is  to  perform  the  reconfiguration 
over  a  user-specified  time  and  meet  the  final  conditions  while  minimizing  the  loss  in  energy 
as  a  result  of  the  maneuver.  Therefore,  for  this  study,  a  linear  quadratic  terminal  control 
law  is  ideal. 

Within  this  chapter,  the  issue  of  controllability  will  be  discussed  for  a  single  satellite 
in  Section  3.1.  Then,  the  state  equations  for  the  controller  will  be  developed  in  Section 
3.2  and  applied  to  the  terminal  controller  in  Section  3.3.  Finally,  the  constraint  equations 
will  be  defined  in  Section  3.4. 
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3.1  State  Equations  and  Controllability 

The  CW  equations  derived  in  Section  2.1.1  will  be  used  to  describe  the  relative 
motion  of  the  satellites  for  this  study.  The  current  definition  relates  the  motion  of  one 
satellite  with  respect  to  a  reference  position,  so  only  one  satellite  will  be  considered  for  the 
moment.  The  only  external  force  acting  on  the  satellite  that  affects  its  relative  motion  with 
respect  to  the  formation  is  atmospheric  drag,  which  acts  against  the  velocity  vector  of  the 
satellite.  Since  this  study  assumes  the  reference  orbit  is  circular,  drag  acts  in  the  negative 
eg  direction,  as  defined  in  Figure  2.1.  It  is  assumed  that  the  relative  velocity  components 
of  the  satellite  do  not  affect  the  drag  acceleration.  Instead,  only  the  orbital  velocity  of  the 
circular  reference  orbit  will  be  considered  since  it  is  on  the  order  of  104  larger  than  the 
relative  velocities.  As  a  result,  the  drag  acceleration  acts  only  on  the  in-track  component 
of  the  motion,  and  all  other  differential  forces  will  be  considered  zero.  Therefore,  the  force 
per  unit  mass  components  from  Equations  2.15-2.17  become 


fr  =  0  (3.1) 

fo  =  ~\B*PV2  (3.2) 

fz  =  0  (3.3) 

The  satellite  is  controlled  by  changing  the  angle  of  attack  of  the  drag  panels  to 
alter  the  presented  area  of  the  satellite.  This  change  in  the  presented  area  is  directly 
proportional  to  the  change  in  the  ballistic  coefficient,  as  shown  in  Equation  2.49.  So 
changing  the  magnitude  of  the  ballistic  coefficient  directly  alters  the  amount  of  drag  acting 
on  the  satellite.  Therefore,  the  ballistic  coefficient  magnitude  is  the  control  input  for  this 
study. 

The  equations  of  motion  for  the  satellite  can  now  be  written  in  state  equation  form 
such  that 


^x(i)  =  Fx(f)  +  Gu(t) 


(3.4) 


3-2 


where 


x 


5r  ?’o 50  5z  5r  r^SO 


(3.5) 


F 


G 


0  0  0 

0  0  0 

0  0  0 

3  n2  0  0 

0  0  0 

0  0  —  n2 


1  0  0 

0  1  0 

0  0  1 

0  2n  0 

—2  n  0  0 

0  0  0 


0  0  0  0 


(3.6) 


(3.7) 


u  =  B* 


(3.8) 


Before  moving  forward  with  these  equations,  the  controllability  should  be  checked. 
For  complete  state  controllability,  an  unconstrained  control  input  that  will  transfer  the 
system  from  some  initial  state,  xo,  to  some  final  state,  xy ,  in  a  finite  time  interval  must 
exist  [16].  From  the  state  equation  definition  shown  in  Equation  3.4,  it  can  be  seen  that 
the  cross-track  component,  5z,  is  decoupled  from  the  in-plane  motion  components,  5r  and 
59.  As  a  result,  any  motion  or  force  that  results  in  the  cross-track  direction  has  no  effect 
on  the  in-plane  motion.  Likewise,  any  motion  or  force  in  the  in-plane  motion  has  no 
effect  on  the  cross-track  motion.  Because  of  this  cross-track  decoupling,  the  5z  component 
cannot  be  controlled  through  drag,  which  acts  on  the  in-plane  motion  only.  This  lack  of 
control  can  be  confirmed  through  the  controllability  matrix.  For  a  system  to  be  completely 
state  controllable,  the  rank  of  the  the  n  x  n  controllability  matrix  must  be  n,  or,  written 
mathematically  [3]: 


Rank 


G  :  FG  :  F2G  : 


Fn~lG  I  =  n 

-  nxn  j 


(3.9) 
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This  condition  verifies  that  each  column  of  the  controllability  matrix  is  linearly  independent 
[16].  If  the  system  is  not  complete  state  controllable,  one  or  more  of  the  states  is  uncon¬ 
trollable.  The  controllability  matrix  for  the  six-state  system  defined  in  Equation  3.4  only 
has  a  rank  of  four,  meaning  two  states  are  uncontrollable.  From  the  decoupled  cross-track 
discussion  above,  we  can  conclude  that  the  dz  and  dz  components  are  indeed  uncontrol¬ 
lable.  Therefore,  redefining  the  state  such  that  only  the  in-plane  motion  is  accounted  for 
yields 


x  = 


dr  ro58  dr  rodO 


(3.10) 


Then,  the  new  system  of  equations  becomes 


dr 

0 

0 

1 

0 

dr 

0 

rod  0 

= 

0 

0 

0 

1 

rod8 

+ 

0 

dr 

3n2 

0 

0 

2  n 

dr 

0 

rod  8 

0 

0 

—2  n 

0 

rod8 

.  -\py2  _ 

(3.11) 


This  new  four-vector  state  yields  a  controllability  matrix  with  a  rank  of  four.  As  a  result, 
the  system  is  completely  controllable. 


3.2  Controller  State  Equations 

Since  the  control  law  is  dependent  on  the  motion  of  both  satellites  in  the  formation, 
the  state  of  the  system  should  be  defined  to  include  the  state  of  each  satellite.  So,  define 
the  controller  state  vector,  x,  as 


such  that 


x(f) 


Xi(f) 

X2  (t) 


8x1 


(3.12) 


x(*)  = 


i  T 


dr  i  r\)d0\  dr  \  rodd\  :  dr 2  rodd-2  dr 2  rod  82 


(3.13) 
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Since  each  satellite  is  controlled  separately,  the  control  input,  u,  is  defined  by 


u{t) 


Bl{t) 


(3.14) 


where  B\  and  are  the  ballistic  coefficients  for  each  satellite  caused  by  the  drag  panel 
orientation.  The  state  equations  then  become 


x(f)  =  Fx(f)  +  Gu{t) 


(3.15) 


where 


Fx 

X 

o 

X 

o 

F2 

G\ 

04x1 

04x1 

g2 

(3.16) 


(3.17) 


Here,  04X4  and  Cbxi  are  (4  x  4)  and  (4x1)  zero  matrices,  respectively.  And  since  the  two 
satellites  are  associated  with  the  same  reference  orbit, 


Fi  =  F2 


0  0  10 

0  0  0  1 

3n2  0  0  2n 

0  0  —2  n  0 


G\  =  G2 


0  0  0 


(3.18) 


(3.19) 


The  matrices  F  and  G  are  assumed  time-invariant  since  they  are  functions  of  the 
orbital  mean  motion,  n,  and  the  satellite  velocity  with  respect  to  the  atmosphere,  v,  of 
the  reference  orbit,  both  of  which  will  be  considered  constant  over  the  controlled  orbit. 
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Because  the  altitude  of  the  reference  orbit  is  decreasing  with  the  controlled  input,  these 
variables  are  actually  changing.  However,  even  for  large  reconfigurations  that  result  in 
a  considerable  change  in  altitude,  the  change  in  mean  motion  and  satellite  velocity  are 
minimal.  As  an  example,  a  large  reconfiguration  with  an  initial  reference  orbital  radius  of 
6800  km  that  causes  an  altitude  loss  of  275  m  results  in  a  0.0061%  change  in  mean  motion 
and  a  0.0025%  change  in  velocity.  Therefore,  the  assumptions  that  F  and  G  are  constant 
are  valid  within  the  bounds  of  this  study. 

The  ballistic  coefficients  considered  for  this  study  are  the  ballistic  coefficients  caused 
by  the  deployment  of  the  drag  panels  only,  not  the  overall  ballistic  coefficients  of  the 
satellites  themselves.  As  discussed  in  Section  2.4,  the  ballistic  coefficients  of  the  satellite 
bodies  will  cancel  in  the  relative  frame,  assuming  the  two  satellites  are  identical.  If  the 
ballistic  coefficients  of  the  satellites  are  not  equal  and  one  satellite  has  a  slightly  larger 
nominal  ballistic  coefficient,  the  second  satellite’s  panels  can  be  oriented  such  that  the  two 
ballistic  coefficients  of  the  satellites  are  equal,  thereby  defining  a  new  nominal  state  for  the 
second  satellite.  As  a  result,  the  second  satellite  will  fly  with  a  higher  drag  acceleration 
than  it  would  normally  sustain.  This  increased  drag  will  shorten  the  life  of  the  formation, 
but  this  tradeoff  is  necessary  to  maintain  the  formation  for  any  extended  period  of  time. 
And  since  only  the  differential  drag  affects  the  motion  of  the  satellites,  the  least  amount  of 
total  drag  should  be  used  to  achieve  the  desired  differential  state.  Therefore,  to  minimize 
the  drag  and,  as  a  direct  result,  the  overall  energy  loss  of  the  formation,  only  one  satellite 
should  have  its  drag  panels  deployed  at  any  time  during  a  reconfiguration.  The  other 
satellite  should  keep  its  panels  in  the  drag  minimizing  state,  which  is  zero  for  the  relative 
frame.  For  example,  when  satellite  1  has  its  drag  panels  deployed,  the  ballistic  coefficient 
of  satellite  2  due  to  its  drag  panels  is  zero  such  that  =  0,  and  vice  versa  when  satellite 
2  has  its  drag  panels  deployed. 

3.3  Linear- Quadratic  Terminal  Controller 

The  goal  of  the  controller  is  to  meet  the  terminal  conditions  while  minimizing  the 
change  in  energy  of  each  satellite  in  the  formation.  As  explained  in  Section  2.3,  since  drag 
acts  against  the  motion  of  the  satellite,  any  change  in  energy  that  results  from  the  control 
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input  in  this  case  decreases  the  overall  energy  of  the  satellite.  Therefore,  minimizing  the 
change  in  energy  is  equivalent  to  minimizing  the  loss  in  energy.  This  characteristic  plays 
well  into  a  quadratic  performance  index  since  squaring  the  change  in  energy  is  equivalent 
to  squaring  the  loss  of  energy.  By  selecting  a  quadratic  performance  index  and  keeping  the 
constraints  linear,  a  linear-quadratic  terminal  controller  can  be  implemented  to  reconfigure 
a  formation  with  optimal  control  inputs  that  are  linear  feedbacks  on  the  state  variables 
[2].  A  performance  index  that  penalizes  the  change  in  orbital  energy  of  each  satellite  is 
chosen,  such  that 


d£  i 

dt 


(3.20) 


where  £\  and  £2  are  the  orbital  energies  of  satellite  one  and  satellite  two,  respectively. 
To  reconfigure  the  formation  to  the  desired  final  formation,  terminal  conditions  must  be 
applied  to  the  control  law.  By  keeping  the  terminal  constraints  soft  such  that  the  terminal 
error  is  driven  close  to  zero  but  not  necessarily  exactly  to  zero,  the  penalty  function  can 
be  made  quadratic  as  well  and  appended  to  the  performance  index.  Define  the  quadratic 
penalty  function  on  the  terminal  error  as 


<P[*(tf)}  =  ^e/Q/e/  (3-21) 

where  Qf  is  an  (8  x  8)  diagonal  error  weighting  matrix,  and  ej  is  the  terminal  error  vector 
given  by 


ef  =  Mfx(tf)  -  if  (3.22) 

Mf  and  if  are  the  terminal  state  matrix  and  the  terminal  state  vector,  respectively,  and 
are  defined  to  satisfy  the  terminal  conditions.  These  two  variables  are  dependent  on  the 
final  formation  and  will  be  defined  in  Section  3.4. 

Appending  the  penalty  function  from  Equation  3.21  to  the  performance  index  from 
Equation  3.20  yields  the  overall  performance  index  [2] 
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(3.23) 


J  ~  2 e/*5/e/  + 


subject  to  the  conditions 


x(t) 

=  .Fx(f)  +  Gu(t) 

(3.24) 

x(i0) 

=  x0 

(3.25) 

Here,  F  and  G  are  the  matrices  defined  in  Equations  3.16  and  3.17  in  the  previous  section, 
and  R  is  the  control  penalty  matrix  that  is  defined  to  satisfy  the  performance  penalty  in 
Equation  3.20.  From  Equations  2.47  and  2.52,  the  energy  rate  of  change  for  one  satellite 
can  be  written  as 


f  -  4b>3  (3-26» 

The  performance  index  defined  in  Equation  3.20  is  quadratic  in  terms  of  the  energy  rate 
of  change.  Therefore,  by  this  definition,  the  control  penalty  matrix  R  must  be  chosen  to 
yield 


d£\ 


dt 


+ 


d£2 


dt 


=  +  I(H2*)VV 


(3.27) 


Since  the  control  input  for  the  two-satellite  system  is  u  =  [HJ,  B%]  ,  R  becomes  the  positive 
definite  matrix  defined  as 


\p2v6  0 

0  \p2v6 


(3.28) 


From  the  performance  index  and  the  equations  of  motion,  the  Hamiltonian,  77,  can 
be  formed,  yielding 
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(3.29) 


TL  =  -uTi?u  +  X1  (fx  +  Gu) 

where  A  is  the  costate  vector.  Then,  from  the  Hamiltonian,  the  Euler-Lagrange  necessary 
conditions  are  formed: 


x 


dH 

~dX 


Fx  +  Gu 


X 


an 

dx 


=  -Ft  X 


0 


an 

du 


Ru  +  GtX 


(3.30) 


(3.31) 


(3.32) 


Solving  Equation  3.32,  which  is  also  known  as  the  stationarity  condition,  for  u (t)  yields 
the  optimal  control  law  in  terms  of  the  costates: 


u  (t)  =  —R~1GrX(t)  (3.33) 

The  hnal  costates  can  be  found  from  the  final  state  penalty  function  defined  in  Equation 
3.21: 


x(tf)  =  =  Mj Qf  [Mfx(tf)  -  ip }  (3.34) 

Then,  combining  Equations  3.30,  3.31,  and  3.33  yields  the  two-point  boundary  value  prob¬ 
lem: 


X 

F  -GR~1Gt 

X 

A 

1 

o 

00 

X 

00 

1 

1 _ 

A 

subject  to  the  boundary  conditions 


(3.35) 
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x(*o)  =  x0 

X(tf)  =  MjQf[Mfx(tf)-il>\ 


(3.36) 

(3.37) 


The  block  diagram  for  this  linear-quadratic  terminal  control  law  is  shown  in  Figure  3.1. 


Figure  3.1  Block  Diagram  of  the  Linear-Quadratic  Terminal  Control  Law 

Since  the  equations  are  linear,  a  state  transition  matrix  can  be  used  to  propagate 
the  final  costates  A/  from  Equation  3.34  backward  in  time  to  obtain  the  initial  costates  Ao 
at  time  to  [2].  Since  the  system  is  also  time- invariant,  the  forward  state  transition  matrix, 
can  be  determined  from  the  matrix  exponential,  such  that  [8] 


Ht,  t0)  =  =  I +  H(t- 10)  +  ^  H\t  -  to)2  +  liH3(t- t0)3  +  •  •  •  +  ^Hk(t  -  t0)k  +  •  •  • 

(3.38) 

where  H  is  the  (16  X  16)  matrix  defined  in  the  two-point  boundary  value  problem  in 
Equation  3.35,  rewritten  here  for  convenience: 


H  = 


F  -GR~1Gt 
08x8  —  FT 


16x16 


The  state  transition  matrix  can  then  be  partitioned  such  that 


(3.39) 
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(3.40) 


x(f) 

$n(L  to) 

$12(Mo) 

x(to) 

A  (t) 

$21  (Mo) 

$22 (t,  to) 

A(i0) 

A  method  of  obtaining  the  initial  costates,  A(to),  now  exists.  From  the  final  costate 
equation  shown  in  Equation  3.37  and  from  the  system  of  equations  defined  in  Equation 
3.40,  A(fo)  becomes 


A(to)  =  {$22 (tf)  -  {  [Mj QfMf$u(tf)  —  $21  (tf)]  x(to)  —  Mj Qf i/>} 

(3.41) 

where  $(t/)  =  4>(t/,to)  for  brevity.  Now  with  the  initial  conditions  x(to)  and  A(fo),  the 
system  of  equations  in  Equation  3.35  can  be  integrated  forward  in  time,  yielding  x(f)  and 
A(f).  Then,  from  Equation  3.33,  the  optimal  control  input  u(f)  can  be  found. 

The  Euler-Lagrange  equations,  shown  in  Equations  3.30-3.32,  are  the  necessary  con¬ 
ditions  that  must  be  satisfied  for  an  optimal  solution.  Meeting  these  conditions  ensures  a 
stationary  solution  has  been  found,  but  it  does  not  necessarily  guarantee  a  minimum  solu¬ 
tion.  A  second-order  sufficient  condition  for  a  minimum  solution  is  given  by  the  Legendre- 
Clebsch  condition,  defined  as  [20] 


5?  > 0  (3-42) 

Satisfying  both  the  Euler-Lagrange  necessary  conditions  and  the  Legendre-Clebsch  suf¬ 
ficient  condition  guarantees  the  control  solution  found  for  u(t)  provides  at  least  a  local 
minimum  for  the  Hamiltonian.  For  this  study,  the  Legendre-Clebsch  condition  yields 

0  =  *  (3.43) 

where  R  is  the  control  penalty  matrix  defined  in  Equation  3.28.  This  matrix  is  positive 
definite,  satisfying  the  Legendre-Clebsch  condition  and  proving  that  the  solution  found  is 
indeed  the  optimal  solution. 
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The  equations  derived  in  this  section  are  dependent  on  the  terminal  state  matrix  Mf 
and  the  terminal  state  vector  -0,  which  dictate  the  final  formation  type.  These  conditions 
will  be  developed  next  in  Section  3.4. 

3-4  Terminal  Constraints 

Since  drag  allows  only  in-plane  control,  the  formation  types  a  drag-actuated  controller 
can  achieve  are  limited.  Therefore,  only  two  formation  designs  will  be  considered  for  this 
study:  in-plane  and  elliptical.  These  formation-types  are  the  bases  for  more  sophisticated 
formations  that  require  cross-track,  or  out-of-plane,  control,  such  as  repeat  groundtrack, 
circular,  and  projected  circular  formations  [18].  The  reconfigurations  of  the  formation 
began  with  some  initial  state  for  the  satellites  based  on  the  initial  formation,  and  the 
terminal  state  for  the  controller  was  defined  for  a  zero  radial  displacement  between  the 
two  such  that,  when  the  control  ended,  the  satellites  had  an  in-track  displacement  only 
with  some  specified  radial  and  in-track  velocity  to  achieve  the  desired  formation.  Also, 
to  eliminate  any  confusion,  the  leading  satellite  in  the  formation  at  the  beginning  of  the 
control  input  has  been  designated  satellite  1,  while  the  trailing  satellite  is  satellite  2.  For 
reconfigurations  where  the  satellites  have  a  zero  in-track  displacement,  the  satellite  at  the 
lower  altitude  is  designated  satellite  1. 

Upon  reconfiguration  of  the  formation,  certain  final  conditions  must  be  satisfied  to 
place  the  satellites  in  the  desired  final  formation.  These  conditions  vary  depending  on  the 
formation  type.  However,  one  condition  that  must  be  met  regardless  of  the  formation  type 
is  a  zero,  or  at  least  a  highly  suppressed,  in-track  drift  between  the  satellites.  Minimizing 
or  zeroing  the  in-track  drift  between  the  two  satellites  keeps  the  uncontrolled  formation 
from  drifting  apart,  reducing  the  need  for  additional  control  to  maintain  the  formation. 
This  reduction  in  control  results  in  long-term  cost  savings,  whether  the  cost  relates  to 
conventional  expendables  or,  in  this  case,  altitude. 

To  suppress  the  in-track  drift,  the  secular  terms  of  the  uncontrolled  solution  devel¬ 
oped  in  Section  2.1.2  must  be  minimized  or  zeroed.  The  in-track  solution  contains  the  only 
secular  terms  in  the  force  free  solution.  From  Equation  2.33,  the  resulting  secular  drift  for 
one  satellite  with  respect  to  the  reference  position  is 
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89 secular  ^  0^^  T  2n8r 


(3.44) 


For  formation  keeping  about  an  empty  reference  position,  the  only  concern  is  that  the  two 
satellites  in  the  formation  are  not  drifting  apart,  not  necessarily  drifting  apart  with  respect 
to  the  reference  point.  So  as  long  as  the  two  satellites  are  drifting  away  from  the  reference 
position  with  the  same  rate  and  direction,  the  formation  will  remain  intact.  Then,  the 
reference  position  can  be  redefined  to  match  the  resulting  drift  from  the  original  reference 
position,  thereby  zeroing  any  drift  in  the  newly  selected  frame.  Therefore,  the  relative  drift 
of  the  satellites  should  be  zeroed,  such  that 

ro(69i  —  862)  +  2n(Sri  —  Sr2 )  =  0  (3.45) 

Achieving  this  constraint  ensures  the  formation  will  not  degrade  due  to  the  nominal  motion 
of  the  satellites. 

Depending  upon  the  desired  formation,  additional  constraints  must  be  imposed  along 
with  the  in-track  constraint  described  above.  The  constraints  required  to  achieve  the  in¬ 
plane  and  elliptical  formations  will  be  discussed. 

3.4.I  In-Plane  Formation.  The  in-plane  formation  is  perhaps  the  simplest  of 
all  formations.  The  two  satellites  share  the  same  orbital  plane  but  are  separated  in  the 
in-track  direction  by  their  true  anomaly  [18].  So  one  satellite  is  flying  ahead  of  the  other 
in  the  same  orbit,  as  shown  in  Figure  3.2.  Again,  since  drag  does  not  allow  cross-track 
control,  the  initial  cross-track  displacement  and  velocity  must  be  zero.  For  the  in-plane 
formation,  a  constraint  must  be  imposed  on  the  radial  component  to  dampen  out  the 
radial  oscillation.  As  the  satellite  radius  increases  with  respect  to  the  reference  position, 
the  satellite’s  orbital  period  increases,  causing  it  to  fall  behind  in  the  in-track  direction. 
Likewise,  when  the  radius  decreases,  the  satellite  moves  ahead  of  the  reference  position. 
As  a  result  of  this  radial  oscillation,  an  in-track  oscillation  also  results.  Therefore,  to 
maintain  the  two  satellites  in  the  same  orbital  plane  and  at  a  fixed  distance  apart,  the 
radial  oscillation  must  be  suppressed.  The  radial  component  from  Equation  2.28  is 
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6r 


Figure  3.2  Illustration  of  the  In-Plane  Formation  in  the  Relative  Frame 


5r(t)  =  —  ( —roSOo  +  3<5?’o  )  cos  nt  +  — ^  sin  nt  +  —roSdo  +  4<5ro  (3.46) 

\n  J  n  n 

The  5tq  and  r^Sd  terms  are  already  being  constrained  to  zero  the  in-track  drift,  as  shown 
in  Equation  3.45.  So  the  only  term  left  to  work  with  is  the  radial  velocity  from  the  sine 
coefficient.  Again,  the  interest  is  to  have  the  two  satellites  move  together.  As  a  result,  the 
sine  coefficients  from  each  satellite  should  be  subtracted,  such  that 

Sri  -Sr2  =  0  (3.47) 

Suppressing  this  term  ensures  the  satellites  will  not  oscillate,  or  will  at  least  oscillate 
together,  in  the  radial  direction  over  time. 

Satisfying  the  constraint  in  Equation  3.47  ensures  the  satellites  are  in  the  same 
orbital  plane,  but  the  distance  between  the  two  satellites  is  not  specified.  Therefore, 
to  place  the  satellites  a  set  distance  apart,  a  constraint  must  be  placed  on  the  in-track 
distance  component  of  the  system  equations.  To  achieve  a  set  separation  distance  of  dsep, 
the  following  constraint  must  be  met: 

roS01  -  r0Sd2  =  dsep  (3.48) 
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This  in-track  separation  distance  constraint,  coupled  with  the  radial  oscillation  constraint 
in  Equation  3.47  and  the  suppressed  in-track  drift  constraint  from  Equation  3.45,  results 
in  the  desired  in-plane  formation. 

These  constraints  can  then  be  used  to  form  the  terminal  constraint  matrix,  M f ,  and 
the  terminal  state  vector,  0,  both  of  which  were  defined  in  Equation  3.22.  Therefore,  from 
the  constraints  defined  in  Equations  3.45,  3.47,  and  3.48,  define  the  very  sparse  diagonal 
matrix  Mf  and  the  terminal  state  vector  0: 


f  (in— plane) 


2  n  0 

010... 

...010 
...  0  1  0 

0  —2  n  0 

0  -1  0  ... 
0-10 
...  0  -1 


(3.49) 


0 


in— plane 


0^00 


0  0 


(3.50) 


Applying  Equations  3.49  and  3.50  in  the  optimal  control  law  defined  in  Section  3.3  results 
in  the  desired  in-plane  formation. 

3-4-2  Elliptical  Formation.  In  an  elliptical  formation,  the  two  satellites  appear  to 
orbit  the  reference  position  on  the  order  of  one  orbital  period  of  the  system.  A  2-by-l  ellipse 
is  the  natural  formation,  with  the  maximum  in-track  separation  being  twice  the  maximum 
radial  separation  as  illustrated  in  Figure  3.3.  These  dimensions  are  evident  in  the  solution 
to  the  CW  equations  shown  in  Equation  2.33.  The  coefficients  of  both  the  cosine  and  the 
sine  terms  of  the  r^SO  solution  are  double  the  coefficients  of  the  Sr  solution,  resulting  in  a 
natural  2-by-l  ellipse  given  any  radial  or  in-track  velocity  or  a  radial  displacement. 
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Figure  3.3  Illustration  of  the  Elliptical  Formation  in  the  Relative  Frame 

For  an  elliptical  formation,  a  prescribed  radial  oscillation  is  required  to  enable  the 
satellite  to  orbit  the  reference  position  in  a  2-by-l  ellipse.  This  prescribed  oscillation  can 
be  achieved  by  just  focusing  on  the  sine  term.  By  choosing  the  proper  radial  velocity 
difference,  the  desired  ellipse  can  be  formed.  Subtracting  the  two  sine  terms  results  in 


- =  Srmax  (3.51) 

n  n 

where  5rmax  is  the  radial  separation  required  to  place  the  satellites  in  the  prescribed 
elliptical  formation.  Since  the  maximum  radial  distance  is  half  the  maximum  in-track 
distance,  Srmax  can  be  solved  for  in  terms  of  Srrei,  the  velocity  magnitude  of  each  satellite 
required  to  place  the  satellites  in  the  elliptical  formation,  such  that 
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.  ndsep 
OTrel  =  ~ ~ 

Then,  the  radial  velocity  constraint  can  be  defined  as 


(3.53) 
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5n  -  Sr2  =  Sfrei 


(3.54) 


By  meeting  this  radial  velocity  constraint  defined  in  Equation  3.54  along  with  the  in-track 
drift  constraint  from  Equation  3.45  and  the  in-track  separation  distance  constraint  from 
Equation  3.48,  each  satellite  will  trace  the  same  2-by-l  ellipse  in  the  radial/in-track  plane. 

These  constraints  can  then  be  used  to  form  the  terminal  state  vector,  -0,  which  was 
defined  in  Equation  3.22.  Since  the  same  states  from  the  in-plane  formation  are  being 
constrained,  the  terminal  constraint  matrix  Mf  is  the  same  as  the  one  defined  in  Equation 
3.49.  Based  on  the  constraints  defined  in  Equations  3.45,  3.48,  and  3.54,  the  terminal 
constraint  matrix  M f  and  the  terminal  state  vector  ■?/>  become 

M f  (elliptical)  Mf  (in— plane)  (3.55) 
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(3.56) 


Applying  this  vector  along  with  the  matrix  defined  in  Equation  3.49  in  the  optimal  control 
law  defined  in  Section  3.3  results  in  the  desired  elliptical  formation. 
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IV.  Results 


Reconfigurations  to  both  the  in-plane  and  elliptical  formations  were  performed  from  initial 
formations  of  in-plane,  elliptical,  and  a  post-payload  separation  position.  The  post-payload 
separation  is  intended  to  represent  an  initial  formation  placement  after  the  satellites  have 
been  released  from  their  launch  vehicle,  assuming  they  were  launched  together.  The  opti¬ 
mal  control  law  presented  in  Chapter  III  is  a  fixed-time  routine,  meaning  the  user  specifies 
the  amount  of  time  over  which  the  optimization  will  take  place.  The  fixed-time  input  allows 
the  user  to  extend  the  time  of  the  reconfiguration  to  reduce  the  control  input  magnitude 
throughout  the  reconfiguration.  The  ability  to  increase  or  decrease  the  time  allows  the 
user  to  specify  a  transfer  that  will  most  effectively  make  use  of  the  allotted  control  input 
without  exceeding  that  amount.  Also,  extending  the  time  allowed  for  the  transfer  tightens 
the  constraints  for  most  cases,  resulting  in  terminal  conditions  that  match  the  terminal 
constraints  more  closely.  A  conservative  value  of  0.1  m2 /kg  was  chosen  for  the  maximum 
differential  ballistic  coefficient  magnitude  for  this  study.  As  for  the  time  allowed  for  the 
reconfiguration,  a  natural  time-scale  for  orbital  motion  is  an  orbital  period.  And  for  the 
elliptical  formation,  each  satellite  completes  one  full  orbit  of  the  reference  position  in  the 
relative  frame  over  one  orbital  period.  As  such,  the  reconfiguration  durations  have  been 
specified  in  terms  of  orbital  periods  for  the  reference  orbit. 

A  concern  when  applying  optimization  routines  with  soft  terminal  constraints  is  the 
final  weight  of  both  the  final  constraints  and  the  control  cost  within  the  performance  index. 
The  goal  is  to  minimize  the  performance  index  while  meeting  the  terminal  constraints 
and  minimizing  the  control  cost.  But  if  too  much  effort  is  applied  to  reducing  the  final 
constraint  error,  the  control  input  is  not  necessarily  optimal.  As  such,  a  small  terminal 
constraint  cost  is  desired,  so  the  final  control  and  terminal  constraint  costs  are  provided 
for  each  maneuver  as  a  performance  measure.  Note  that  the  units  are  ignored  for  these 
costs  since  only  the  weight  of  each  measure  is  needed.  Also,  when  dealing  with  spacecraft 
transfers  and  formation  reconfigurations,  maneuver  costs  are  generally  quantified  by  a  Av 
budget  since  the  actual  fuel  budget  is  dependent  upon  the  spacecraft  mass.  While  the 
method  of  transfer  for  this  study  is  different  from  conventional  techniques,  a  Av  budget 
will  be  determined  for  each  reconfiguration  based  on  the  total  amount  of  control  applied 
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to  each  satellite.  Due  to  the  continuous,  low-thrust  nature  of  the  control  input,  these 
maneuvers  will  result  in  a  higher  fuel  cost  than  optimal  impulsive  burn  maneuvers  using 
conventional  thrusters,  but  the  costs  are  within  the  same  order  of  magnitude.  Also,  these 
costs  are  associated  with  different  fuel  measures;  a  conventional  thruster  fuel  cost  is  related 
directly  to  expendables  onboard  the  spacecraft,  while  the  fuel  cost  in  this  study  is  related 
to  altitude.  Despite  the  difference,  a  Av  budget  is  provided  for  each  reconfiguration. 

The  post  reconfiguration  plots  presented  in  Sections  4.2  and  4.3  contain  a  radial 
position  deviation  from  the  expected  final  position  and  an  in-track  separation  between  the 
two  satellites  over  the  duration  of  the  20  orbital  period  propagation.  The  radial  position 
deviation  is  on  the  sub- millimeter  level  for  most  cases,  a  level  of  precision  that  is  most  likely 
undetectable  with  current  positioning  techniques.  But  the  author  feels  it  is  important  to 
present  this  deviation  to  provide  insight  into  the  cause  of  the  in-track  drift,  which  is  on  the 
centimeter  level  over  the  duration  of  the  propagation.  This  centimeter-level  drift  would 
be  detectable  for  phased  array  applications,  so  for  a  thorough  investigation,  the  drift  is 
presented. 

This  chapter  presents  the  results  of  the  simulation  based  on  the  equations  of  motion 
presented  in  Chapter  II  and  the  control  law  presented  in  Chapter  III.  Section  4.1  verifies 
the  equations  of  motion  used  for  this  study  were  properly  coded  by  checking  against  a 
separate  set  of  equations.  Sections  4.2  and  4.3  present  the  results  of  reconfigurations  into 
in-plane  and  elliptical  formations.  All  of  the  results  presented  within  this  chapter  had 
a  reference  orbital  radius  of  6800  km.  For  the  reverse  of  the  reconfigurations  presented 
within  this  chapter,  see  Appendix  A,  and  for  reconfigurations  performed  at  a  reference 
orbital  radius  of  7000  km,  see  Appendix  B. 

4-1  Code  Verification 

The  equations  of  motion  were  numerically  integrated  using  Matlab’s  ode45  inte¬ 
grator.  To  verify  the  integrated  equations  of  motion  were  correct,  both  the  unforced  and 
the  forced  equations  were  tested  using  different  methods.  The  unforced  equations  were 
verified  versus  the  force  free  solution  presented  in  Section  2.1.2.  Integrating  over  the  same 
final  time  as  that  entered  for  the  final  time  of  the  closed  form  solution  given  a  zero  control 
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input  and  equal  initial  conditions,  the  solutions  should  be  equal  within  integration  error. 
For  the  forced  solution,  a  state  transition  matrix  solution  was  formed  to  verify  the  results 
from  the  numerical  integration.  Both  of  these  validation  tests  will  be  further  discussed. 

4-1.1  Unforced  Solution.  To  verify  the  numerically  integrated  solution  to  the 
unforced  equations  of  motion,  the  closed  form  solution  developed  in  Section  2.1.2  was  used. 
The  equations  of  motion  of  a  satellite  was  integrated  over  three  orbital  periods.  With  a 
reference  orbital  radius  of  6800  km,  the  state  vector  of  the  satellite  for  this  example  was 


xq 


T 


10  m  0  0  0 


(4.1) 


The  total  time  was  divided  into  500  equal  intervals,  and  the  differences  between  the  radial 
and  in-track  results  for  both  the  numerically  integrated  solution  and  the  closed  form  so¬ 
lution  were  plotted.  Figure  4.1  shows  the  plot  of  the  solution  in  the  radial/in-track  plane 
to  give  an  idea  of  the  general  motion  of  the  satellite,  and  Figure  4.2  shows  the  radial  and 
in-track  differences  over  the  length  of  the  integration.  As  shown  in  Figure  4.1,  the  satellite 
falls  behind  the  reference  position  as  expected  due  to  its  larger  orbital  radius,  which  results 
in  a  larger  orbital  period.  The  radius  also  oscillates  as  a  direct  result  of  its  slight  orbital 
eccentricity,  returning  to  its  original  position  with  each  orbital  period.  Figure  4.2  shows 
the  two  solutions  agree  to  the  sub-micrometer  level  over  the  three  orbital  periods,  proving 
the  equations  of  motion  are  properly  coded.  The  solutions  begin  to  deviate  more  toward 
the  end  of  the  integration  period  due  to  the  integration  error,  not  errors  in  the  equations. 
Because  the  integrator  uses  only  the  previous  solution  for  each  subsequent  time-step,  dis¬ 
carding  all  other  points,  the  integration  error  tends  to  build  with  time.  However,  the  error 
is  acceptable  for  this  study  since  the  maximum  time  of  integration  in  most  cases  is  limited 
to  20  orbital  periods,  although  this  time  is  increased  for  several  simulations  with  an  orbital 
radius  of  7000  km. 


4-1.2  Forced  Solution.  While  the  force  free  solution  to  the  CW  equations  has  a 
nice  closed  form  solution,  such  an  analytical  solution  does  not  exist  for  the  forced  equa¬ 
tions.  Even  the  forced  solution  presented  in  Section  2.1.3  must  be  integrated.  Also,  the 
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Figure  4.1  In-Plane  Plot  of  the  Unforced  Solution 
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Figure  4.2  Radial  and  In- Track  Differences:  Unforced  Solution 


way  in  which  the  forced  solution  is  presented  for  this  study  is  slightly  different  from  the 
solution  shown  in  Section  2.1.3.  The  initial  control  input  in  terms  of  the  costates,  X(to),  is 
prescribed,  and  the  following  control  inputs  are  then  determined  from  the  overall  system 
of  equations.  To  correctly  verify  the  numerically  integrated  system  of  equations,  a  system 
that  performs  in  the  same  way  should  be  implemented.  Therefore,  a  separate  linear  model 
was  developed  to  meet  this  requirement. 

Using  the  same  state  transition  matrix  derived  in  Section  3.3,  an  approximate  linear 
solution  is  developed,  such  that 


x(£  +  At) 

(i  +  At,  t)  $i2(t  +  At,t) 

x(t) 

A  (t  +  At) 

$21  (t  +  At ,  t)  $22  (t  +  At ,  t) 

.  A(i)  _ 

As  can  be  seen  from  this  solution,  the  subsequent  inputs  for  both  the  state,  x(t),  and 
costates,  A(t),  are  contained  within  the  system,  equivalent  to  the  system  of  equations  that 
are  integrated.  As  such,  this  linear  solution  should  provide  a  solution  that  relates  closely 
to  the  numerical  integration. 


In-Track  (m) 


Figure  4.3  In-Plane  Plot  of  the  Forced  Solution 
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Figure  4.4  Radial,  In- Track,  and  Control  Input  Differences:  Forced  Solution 

To  check  the  numerically  integrated  system  of  equations  against  this  linear  formula¬ 
tion,  a  reconfiguration  was  ordered  from  an  initial  in-track  formation  with  a  50  m  separation 
to  a  final  elliptical  formation  with  an  in-track  width  of  200  m.  The  entire  transfer  will  not 
be  performed;  instead,  only  the  first  three  orbital  periods  will  be  calculated  to  remain 
consistent  with  the  unforced  solution  verification.  Propagating  over  three  orbital  periods 
provides  a  sufficient  number  of  data  points  to  verify  the  equations  are  being  implemented 
correctly.  Figure  4.3  shows  the  in-plane  motion  of  the  two  satellites  in  the  falling  relative 
frame  along  with  their  control  inputs  with  respect  to  time.  These  plots  were  provided  to 
give  an  idea  of  the  overall  motion  and  inputs  of  each  satellite  in  the  formation,  with  the 
motion  of  Satellite  1  depicted  by  the  solid  line  and  Satellite  2  by  the  dotted  line.  Figure  4.4 
shows  the  differences  between  the  numerically  integrated  solution  and  the  linear  solution 
for  the  radial  displacement  and  in-track  displacement,  and  plots  the  percent  difference  for 
the  control  input  solutions.  As  shown  by  the  plots,  the  radial  solutions  agree  to  within  a 
few  millimeters,  while  the  in-track  solutions  agree  to  sub-millimeter  accuracy.  Also,  the 
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control  inputs  agree  on  the  order  of  10  6  of  a  percent,  verifying  that  the  forced  system  of 
equations  are  properly  coded. 

4-2  Reconfiguration  to  an  In- Plane  Formation 

This  section  contains  the  results  for  reconfigurations  into  a  final  in-plane  formation 
with  a  separation  distance  of  500  m  from  three  separate  initial  formations:  an  in-plane 
formation  with  a  separation  distance  of  50  m,  an  elliptical  formation  with  a  semimajor 
axis  of  25  m  such  that  the  satellites  are  50  m  apart  in  the  in-track  position,  and  a  post¬ 
payload  separation  with  initial  conditions  defined  to  represent  two  satellites  that  have 
recently  separated  upon  delivery  to  their  parking  orbit.  The  satellite  positions  and  motion 
are  with  respect  to  a  6800  km .  reference  orbit. 

4-2.1  In-Plane  to  In- Plane.  An  initial  in-plane  formation  with  a  separation 

distance  of  50m.  was  reconfigured  to  a  final  in-plane  formation  with  a  separation  distance 
of  500  m.  This  reconfiguration  is  relatively  easy  to  perform  since  both  satellites  currently 
have  no  radial  oscillations.  So  the  goal  for  the  reconfiguration  is  to  induce  an  in-track  drift 
to  increase  the  size  of  the  formation,  then  zero  the  drift  by  bringing  the  satellites  back 
in-line  at  the  same  altitude.  Since  this  is  a  relatively  simple  reconfiguration,  only  four 
orbital  periods  were  required  for  the  transfer.  Figure  4.5  shows  the  control  inputs  for  each 
satellite,  with  satellite  1  plotted  with  a  solid  line  and  satellite  2  plotted  with  a  dotted  line. 
A  maximum  ballistic  coefficient  input  of  only  0.029  m2/kg  was  required  to  perform  the 
transfer,  well  below  the  maximum  value  allowed.  To  effect  the  reconfiguration,  satellite 
1,  the  leading  satellite,  deployed  its  panels  first,  dropping  in  altitude  to  gain  the  needed 
separation  before  satellite  2  initiated  its  control  input.  The  figure  also  shows  that  both 
satellites  required  the  same  control  input,  just  in  reverse  order,  resulting  in  a  final  Av 
of  0.01048  m/sec  for  each  satellite  to  perform  the  maneuver.  Figures  4.6  and  4.7  show 
the  radial/in-track  position  plots  of  the  satellite  motion  in  both  the  fixed  relative  frame 
and  the  falling  relative  frame.  As  shown  in  the  fixed  frame,  the  overall  altitude  decreases 
approximately  18  m  due  to  the  reconfiguration.  The  plots  also  confirm  the  control  usage 
as  both  satellites  follow  equivalent  trajectories  to  their  final  states.  Figure  4.8  shows  the 
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radial  and  in-track  displacements  from  the  reference  position  plotted  over  the  duration 
of  the  reconfiguration.  This  figure  shows  that  the  satellites  exhibit  acceptable  behavior 
throughout  the  transfer  and  that  the  motion  of  the  satellites  responded  well  to  the  input. 
Also,  the  maximum  radial  displacement  occurs  midway  through  the  reconfiguration,  which 
is  expected  to  occur  at  this  point  since  the  two  satellites  increase  their  relative  velocities 
over  the  first  half  of  the  reconfiguration,  and  then  use  the  second  half  of  the  transfer  to 
decrease  their  relative  velocities  to  zero  as  they  approach  the  500  m  separation  distance. 
Finally,  the  optimality  performance  was  verified,  with  the  performance  index  yielding  a 
control  cost  of  0.65,  while  the  terminal  constraint  cost  was  only  le  —  5,  several  orders  of 
magnitude  less. 

After  the  transfer  was  completed,  the  post-reconfiguration  analysis  was  performed  at 
the  new  altitude  and  with  the  new  mean  motion.  The  final  state  from  the  controlled  motion 
was  propagated  over  20  orbital  periods,  or  approximately  30  hours,  to  ensure  the  transfer 
met  the  desired,  final  formation.  Figure  4.9  shows  the  result  of  the  reconfiguration  in  the 
radial/in-track  plane.  As  shown  in  the  figure,  a  final  in-plane  formation  with  an  in-track 
separation  distance  of  500  m  resulted  from  the  transfer.  The  figure  contains  the  radial 
and  in-track  position  for  each  satellite  over  the  full  20  orbital  periods.  From  the  figure, 
no  drift  is  evident  between  the  two  satellites,  but  due  to  the  scale,  a  small  drift  would 
be  difficult  to  detect.  Therefore,  the  difference  between  the  actual  and  expected  radial 
position,  as  well  as  the  difference  between  the  actual  and  expected  in-track  separation 
between  the  two  satellites,  were  calculated  for  this  position  data  set  and  were  plotted  over 
the  duration  over  the  propagation,  as  shown  in  Figure  4.10.  As  evident  from  the  figure,  each 
satellite  has  a  slight  oscillation  in  the  radial  direction  with  an  amplitude  of  0.24  mm,  and 
the  two  satellites  have  slightly  different  orbital  radii,  with  satellite  2  flying  approximately 
0.5  mm  above  satellite  1  on  average.  This  radial  difference  leads  to  a  slight  drift  between 
the  two  satellites,  causing  the  formation  to  spread  approximately  8  cm  over  the  20  orbit 
propagation,  as  shown  in  the  in-track  deviation  plot.  This  drift  is  well  within  the  bounds 
of  the  constraints  and  does  not  have  noticeable  effects  on  the  formation,  proving  that  the 
optimization  routine  satisfied  the  constraints  of  the  problem. 
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Radial  (m)  Control  Input:  BC  (m2/kg) 


Table  4.1  In-Plane  to  In-Plane  Reconfiguration  Details,  6800  km 


Reference  Orbital  Radius 

6800  km 

Initial  Formation 

In-Plane:  50  m  separation 

Final  Formation 

In-Plane:  500  m  separation 

Total  Time  for 

4  Orbital  Periods 

Reconfiguration 

(6  hr,  12  min,  2  sec ) 

Maneuver 

Avi  =  0.01048  ml  sec 

Cost 

Av2  =  0.01048  ml  sec 

Altitude  Loss 

18.41  m 

Control  Cost 

0.65 

Constraint  Cost 

le— 5 

Figure  4.5  Control  Input:  In-Plane  to  In-Plane,  6800  km 


Figure  4.6  Fixed  Frame  Plot:  In-Plane  to  In-Plane,  6800  km 
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In-track  (m)  Radial  (m)  Radial  (m) 


Figure  4.7  Falling  Frame  Plot:  In-Plane  to  In-Plane,  6800  km 
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Figure  4.8  Reconfiguration  Time  History  Plot:  In-Plane  to  In-Plane,  6800  km 


4-2.2  Elliptical  to  In-Plane.  The  next  reconfiguration  performed  started  from  an 
elliptical  formation  with  an  in-track  separation  distance  of  50  m,  or  a  relative  semimajor 
axis  of  25  m,  and  transferred  to  an  in-plane  formation  with  a  separation  distance  of  500  m, 
the  same  final  formation  as  in  Section  4.2.1.  To  effect  the  transfer,  the  controller  must 
reduce  the  radial  oscillations  of  both  satellites,  effectively  zeroing  their  eccentricities,  while 
driving  them  into  the  proper  in-plane  formation  at  the  correct  separation  distance.  Due 
to  the  initial  orbital  eccentricities  of  the  satellites,  this  transfer  required  a  slightly  larger 
control  input  than  the  previous  reconfiguration.  So  the  controlled  time  was  increased  to  six 
orbital  periods  for  this  maneuver.  Figure  4.11  shows  the  control  inputs  for  each  satellite. 
With  the  increased  time  for  reconfiguration,  a  maximum  ballistic  coefficient  input  of  only 
0.028  m2 /kg  was  required  for  the  transfer.  Again,  notice  that  the  control  inputs  for  the  two 
satellites  are  equal,  but  in  reverse  order.  This  trend  is  expected  since  the  satellites  have  the 
same  initial  conditions  at  the  start  of  the  reconfiguration,  but  the  directions  are  reversed. 
As  a  result,  both  satellites  required  the  same  total  Av  of  0.01151  m/sec  to  perform  the 
reconfiguration.  Figures  4.12  and  4.13  show  the  reconfiguration  in  both  the  fixed  and 
falling  reference  frames.  As  shown  in  the  fixed  frame,  the  satellites  drop  19  m  in  altitude 
as  a  result  of  the  reconfiguration.  These  plots  also  give  insight  into  the  trajectory  taken 
by  each  satellite.  Both  satellites  dampened  their  radial  oscillation  as  they  approached 
the  appropriate  in-track  separation  distance,  but  satellite  1  initially  decreased  its  radial 
position  while  satellite  2  increased  its  radius  to  effect  the  drift,  taking  full  advantage  of 
their  natural  dynamics.  Figure  4.14  confirms  this  approach,  showing  the  time  history  plots 
of  the  reconfiguration  for  both  the  radial  and  the  in-track  distance  components.  Again, 
both  satellites  exhibit  acceptable  behavior  throughout  the  reconfiguration.  The  radial 
oscillation  is  dampened  to  zero  as  the  in-track  separation  distance  is  increased  500  m.  The 
transfer  resulted  in  a  control  cost  of  0.61,  while  the  terminal  constraint  cost  was  only  4e  — 6, 
again  several  orders  of  magnitude  less. 

The  post-reconfiguration  results  were  then  propagated  over  20  orbital  periods  to  show 
the  resulting  motion  of  the  formation.  Figure  4.15  shows  a  plot  of  the  satellites’  positions  in 
the  radial/in-track  plane  over  the  duration  of  the  propagation.  The  satellites  maintain  the 
in-plane  formation  with  a  separation  distance  of  500  m.  without  much  deviation.  Plots  of 
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the  radial  and  in-track  position  deviations  over  time  are  shown  in  Figure  4.16.  Again,  the 
radial  velocity  components  were  not  driven  completely  to  zero,  resulting  in  sub- millimeter 
oscillations  with  each  orbital  period,  well  within  the  bounds  of  the  constraints  but  still 
enough  to  result  in  a  slight  drift.  Since  satellite  2  has  a  larger  average  radial  displacement, 
a  positive  drift  in  the  formation  results,  causing  the  satellites  to  drift  apart  4.5  cm  over  the 
20  orbital  periods. 

Table  4.2  Elliptical  to  In-Plane  Reconfiguration  Details,  6800  km 


Reference  Orbital  Radius 

6800  km 

Initial  Formation 

Elliptical:  25  m  semimajor  axis 

Final  Formation 

In-Plane:  500  m  separation 

Total  Time  for 

6  Orbital  Periods 

Reconfiguration 

(9  hr ,  18  min ,  3  sec) 

Maneuver 

Avi  =  0.01151  m/sec 

Cost 

Av2  =  0.01151  ml  sec 

Altitude  Loss 

20.12  m 

Control  Cost 

0.61 

Constraint  Cost 

4e— 6 

Figure  4.11  Control  Input:  Elliptical  to  In-Plane,  6800  kin 
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Figure  4.12  Fixed  Frame  Plot:  Elliptical  to  In-Plane,  6800  km 


Figure  4.13  Falling  Frame  Plot:  Elliptical  to  In-Plane,  6800  km 


Figure  4.14  Reconfiguration  Time  History  Plot:  Elliptical  to  In-Plane,  6800  fern 
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Figure  4.15  Post  Transfer  Plot:  Elliptical  to  In-Plane,  6800  km. 


Figure  4.16  Post  Transfer  Time  History:  Elliptical  to  In-Plane,  6800  km 


4-2.3  Post-Payload  Separation  to  In-Plane.  Finally,  a  transfer  to  an  in-plane 
formation  from  a  post-payload  separation  configuration  was  performed.  The  initial  states 
for  the  two  satellites  were  defined  as 


xi  =  [—10m,  20m,  — 0.04m/s,  0.005m/s]  (4.3) 

X2  =  [10m,  —20m,  0.04m/s,  — 0.005m/s]  (4.4) 

Again,  the  final  in-plane  formation  had  a  separation  distance  of  500  m,  equal  to  the 
previous  two  reconfigurations  as  shown  in  Table  4.3.  The  results  from  this  transfer  are  quite 
different  from  the  previous  two  reconfigurations  to  an  in-plane  formation.  First,  from  the 
control  input  shown  in  Figure  4.17,  satellite  2  uses  more  control  than  satellite  1,  whereas 
the  previous  transfers  required  equal  control  inputs  for  each  satellite.  Some  intuition  into 
this  control  difference  can  be  seen  in  the  plot  of  the  fixed  reference  frame,  shown  in  Figure 
4.18.  Satellite  1  takes  advantage  of  its  natural  dynamics  from  the  initial  conditions,  going 
uncontrolled  for  most  of  three  orbital  periods,  allowing  its  unforced  trajectory  to  carry  it 
forward  approximately  1200  m  before  slowing  itself  down  and  dampening  its  oscillation. 
The  natural  dynamics  of  satellite  2,  on  the  other  hand,  were  carrying  it  in  the  opposite 
direction,  requiring  the  satellite  to  expend  substantial  control  early  into  the  reconfiguration 
to  overcome  this  nominal  motion.  As  a  result,  satellite  2  required  a  fuel  budget  of  At>2  = 
0.05644  m/sec,  while  satellite  1  only  required  a  budget  of  Aiq  =  0.02139  m/sec,  less  than 
half  the  amount  of  satellite  2.  Despite  the  large  reconfiguration,  the  maximum  ballistic 
coefficient  input  required  was  only  0.06  w? /kg  due  to  the  amount  of  time  allowed  for  the 
reconfiguration.  Figure  4.20  shows  the  time  history  plots  of  the  radial  and  in-track  positions 
over  the  duration  of  the  reconfiguration.  Again,  both  satellites  exhibit  acceptable  behavior, 
dampening  out  their  radial  oscillations  over  the  course  of  the  reconfiguration.  Because  of 
the  larger  reconfiguration  and  the  need  for  a  larger  control  input,  the  performance  index 
cost  was  much  higher,  yielding  a  control  cost  of  4.14,  while  the  constraint  cost  was  lower 
than  the  previous  maneuvers  at  only  2e  —  6. 
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The  post-reconfiguration  propagation  results  can  be  seen  in  Figures  4.21  and  4.22. 
The  radial/in-track  position  plot  shows  the  satellites  in  the  in-plane  formation  with  a 
separation  distance  of  500  m  without  much  deviation.  The  position  deviation  plots  in 
Figure  4.22  shows  a  small  radial  oscillation  with  a  total  displacement  of  approximately 
1.5  mm  that  is  slightly  offset  from  zero,  but  both  satellites  are  nearly  oscillating  together 
with  satellite  1  following  a  slightly  larger  orbit.  As  a  result,  a  negative  drift  results,  causing 
the  satellites  to  move  together  2  cm  over  the  propagation.  Again,  these  results  are  well 
within  the  bounds  of  the  constraints  and  show  that  the  optimization  routine  met  the 
desired  results. 

Table  4.3  Payload  Separation  to  In-Plane  Reconfiguration  Details,  6800  km 


Reference  Orbital  Radius 

6800  km 

Initial  Conditions 

xi  =  [—10m,  20m,  — 0.04m/s,  0.005  m/s] 
X2  =  [10  m,  —20  m,  0.04  m/s,  —  0.005  m/s] 

Final  Formation 

In-Plane:  500  m  separation 

Total  Time  for 

10  Orbital  Periods 

Reconfiguration 

(15  hr,  30  min,  5  sec) 

Maneuver 

Avi  =  0.02139  ml  sec 

Cost 

Av2  =  0.05644  ml  sec 

Altitude  Loss 

64.65  m 

Control  Cost 

4.14 

Constraint  Cost 

2e— 6 

Figure  4.17  Control  Input:  Payload  Separation  to  In-Plane,  6800  km 
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Figure  4.18  Fixed  Frame  Plot:  Payload  Separation  to  In-Plane,  6800  km 
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Figure  4.19  Falling  Frame  Plot:  Payload  Separation  to  In-Plane,  6800  km 
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Figure  4.20  Reconfiguration  Time  History  Plot:  Payload  Separation  to  In-Plane, 
6800  km 
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Figure  4.21  Post  Transfer  Plot:  Payload  Separation  to  In-Plane,  6800  km 
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Figure  4.22  Post  Transfer  Time  History:  Payload  Separation  to  In-Plane,  6800  km 


4-3  Reconfiguration  to  an  Elliptical  Formation 

This  section  contains  the  results  for  reconfigurations  into  a  final  elliptical  formation 
with  a  semimajor  axis  of  250  m,  such  that  the  satellites  are  500  m  apart  when  aligned  in 
the  in-track  direction,  from  three  separate  initial  formations:  an  in-plane  formation  with 
a  separation  distance  of  50  m,  an  elliptical  formation  with  a  semimajor  axis  of  25  m,  and 
a  post-payload  separation  with  the  same  initial  conditions  defined  in  Section  4.2.  The 
satellite  positions  and  motion  are  with  respect  to  a  6800  km  reference  orbit. 

4-3.1  In-Plane  to  Elliptical.  The  first  transfer  performed  reconfigured  a  two- 
satellite  in-plane  formation  with  a  separation  distance  of  50  m  to  an  elliptical  formation 
with  a  semimajor  axis  of  250  m.  The  goal  here  is  to  create  a  radial  oscillation  that  will 
project  the  satellites  into  the  desired  formation  while  limiting  the  drift  and  defining  the 
proper  dimensions  of  the  relative  orbit.  Creating  the  radial  oscillation  to  project  the 
satellites  into  the  correct  orbital  eccentricity  required  more  control  effort  than  the  previous 
reconfigurations,  so  the  transfer  was  performed  over  ten  orbital  periods  to  reduce  the 
maximum  input  required.  Figure  4.23  shows  the  control  input  required  to  effect  the  transfer 
for  both  satellites.  A  maximum  ballistic  coefficient  input  of  0.084  m2 /kg  was  required  for 
each  satellite.  Again,  the  inputs  are  equivalent,  but  in  reverse  order  to  most  effectively 
take  advantage  of  the  natural  motion  of  the  system.  This  total  control  input  resulted  in  a 
Av  for  each  satellite  of  0.08964  m/sec.  The  fixed  and  falling  reference  frame  plots  of  the 
radial/in-track  motion  of  the  satellites  can  be  seen  in  Figures  4.24  and  4.25.  As  shown  in 
the  fixed  frame  plot,  the  formation  altitude  decreases  approximately  150  m  as  a  result  of 
the  reconfiguration.  The  reference  frame  plots  show  that  the  satellites  take  advantage  of 
their  natural  motion,  gradually  increasing  their  radial  oscillation  as  they  spiral  outward 
to  the  desired  position.  To  gain  the  needed  separation,  satellite  1  decreased  its  average, 
relative  altitude  as  satellite  2  increased  its  average,  relative  altitude  to  obtain  the  needed 
in-track  separation.  These  dynamics  are  also  evident  in  the  radial  and  in-track  position 
time  history  plots  shown  in  Figure  4.26.  The  radial  oscillation  increases  over  the  duration 
of  the  reconfiguration,  producing  the  needed  in-track  oscillation  as  well.  The  plots  show 
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that  the  satellites  exhibit  acceptable  behavior  throughout  the  transfer,  remaining  within 
the  bounds  of  the  final  motion. 

The  resulting  state  from  the  reconfiguration  was  then  propagated  over  20  orbital 
periods  to  ensure  the  desired  conditions  were  met.  The  radial/in-track  position  plot  of 
the  satellite  motion  is  shown  in  Figure  4.27.  As  shown  in  the  figure,  the  desired  elliptical 
formation  resulted  from  the  reconfiguration  with  both  satellites  following  the  same  relative 
orbit  with  no  noticeable  drift.  The  position  deviation  plots  in  Figure  4.28  shows  the  radial 
oscillation  is  slightly  larger  than  prescribed,  but  well  within  any  reasonable  constraint 
bounds.  Also,  satellite  1  is  in  a  slightly  larger  orbit,  effecting  a  2  cm  negative  drift  between 
the  satellites  over  the  propagated  time. 


Table  4.4  In-Plane  to  Elliptical  Reconfiguration  Details,  6800  /cm. 


Reference  Orbital  Radius 

6800  km 

Initial  Formation 

In-Plane:  50  m  separation 

Final  Formation 

Elliptical:  250  m.  semimajor  axis 

Total  Time  for 

10  Orbital  Periods 

Reconfiguration 

(15  hr,  30  min,  5  sec) 

Maneuver 

Ai>i  =  0.08964  ml  sec 

Cost 

Av2  =  0.08964  ml  sec 

Altitude  Loss 

148.92  m 

Control  Cost 

18.32 

Constraint  Cost 

le-5 

Figure  4.23  Control  Input:  In-Plane  to  Elliptical,  6800  km 
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Figure  4.24  Fixed  Frame  Plot:  In-Plane  to  Elliptical,  6800  km 


figure  4.25  Falling  Frame  Plot:  In-Plane  to  Elliptical,  6800  km 


Figure  4.26  Reconfiguration  Time  History  Plot:  In-Plane  to  Elliptical,  6800  km 
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Figure  4.27  Post  Transfer  Plot:  In-Plane  to  Elliptical,  6800  km. 


Figure  4.28  Post  Transfer  Time  History:  In-Plane  to  Elliptical,  6800  km 


4-3.2  Elliptical  to  Elliptical.  An  elliptical  formation  with  a  semimajor  axis  of 
25  to  was  expanded  to  250  m  for  this  reconfiguration.  Both  satellites  already  have  a  radial 
oscillation,  so  the  magnitude  of  the  oscillation  must  be  increased  while  the  satellites  move 
apart  to  position  themselves  into  the  new  formation.  This  reconfiguration  was  performed 
over  8  orbital  periods,  as  shown  in  Table  4.5,  which  provides  an  overview  of  the  transfer. 
The  control  inputs  for  each  satellite  are  shown  in  Figure  4.29.  Again,  both  control  inputs 
are  equal,  resulting  in  a  total  Av  of  0.08069m/sec  for  each  satellite  with  a  maximum 
ballistic  coefficient  input  of  approximately  0.096 m? /kg.  Figures  4.30  and  4.31  show  the 
radial/in-track  position  plots  of  the  satellite  motion  in  both  the  fixed  and  falling  frames. 
From  the  falling  frame,  it  is  evident  that  the  controller  is  simply  spiralling  the  two  satellites 
outward,  increasing  both  their  radial  oscillation  and  in-track  separation  distance  over  the 
length  of  the  reconfiguration.  Figure  4.32  confirms  this  motion,  showing  a  radial  and  in¬ 
track  oscillation  whose  amplitude  increases  with  time.  The  radial  and  in-track  position  time 
histories  are  very  similar  for  the  two  satellites,  both  increasing  to  the  final  displacements 
over  the  controlled  time  without  exceeding  the  final  bounds.  Staying  within  these  bounds 
confirms  that  no  control  was  wasted  as  the  controller  allocated  just  enough  to  achieve 
the  desired  result.  The  performance  index  yielded  a  control  cost  of  18.6  with  a  terminal 
constraint  cost  of  2e  —  5. 

Again,  the  final  state  from  the  reconfiguration  was  propagated  over  20  orbital  periods 
to  ensure  the  final  conditions  were  satisfied.  Figure  4.33  shows  the  plot  of  the  radial/in¬ 
track  motion  over  the  propagated  time.  As  evident  from  the  plot,  the  resultant  elliptical 
formation  has  a  semimajor  axis  of  250m,  as  desired,  without  any  noticeable  drift.  More¬ 
over,  the  two  satellites  are  tracing  the  same  relative  trajectory  without  any  deviation.  The 
position  deviation  plots  shown  in  Figure  4.34  further  confirm  that  the  two  satellites  are 
in  the  proper  formations,  although  differing  from  the  prescribed  on  the  millimeter  level. 
The  resulting  semimajor  axis  of  satellite  1  is  smaller  than  that  of  satellite  2  on  the  sub- 
millimeter  level,  effecting  a  negative  drift  of  3  cm  over  the  course  of  the  propagation.  These 
differences  are  again  well  within  the  bounds  of  the  constraints,  proving  the  reconfiguration 
was  successful. 
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Table  4.5  Elliptical  to  Elliptical  Reconfiguration  Details,  6800  km 


Reference  Orbital  Radius 

6800  km 

Initial  Formation 

Elliptical:  25  m  semimajor  axis 

Final  Formation 

Elliptical:  250  m  semimajor  axis 

Total  Time  for 

8  Orbital  Periods 

Reconfiguration 

(12  hr,  24  mm,  4  sec) 

Maneuver 

Avi  =  0.08069  ml  sec 

Cost 

Av2  =  0.08069  ml  sec 

Altitude  Loss 

134.05  m 

Control  Cost 

18.60 

Constraint  Cost 

2e— 5 

Figure  4.29  Control  Input:  Elliptical  to  Elliptical,  6800  Atm 


Figure  4.30  Fixed  Frame  Plot:  Elliptical  to  Elliptical,  6800  km 
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Figure  4.33  Post  Transfer  Plot:  Elliptical  to  Elliptical,  6800  km 


Figure  4.34  Post  Transfer  Time  History:  Elliptical  to  Elliptical,  6800  km 


4-3.3  Post-Payload  Separation  to  Elliptical.  Finally,  a  transfer  to  an  elliptical 
formation  from  a  post-payload  separation  configuration  was  performed.  The  initial  states 
for  this  reconfiguration  were  the  same  used  for  the  reconfiguration  to  an  in-plane  formation, 
but  are  reprinted  here  for  convenience: 


xi  =  [—10m,  20m,  —  0.04m/s,  0.005m/s]  (4.5) 

X2  =  [10m,  —20m,  0.04m/s,  — 0.005m/s]  (4.6) 

The  final  elliptical  formation  had  a  semimajor  axis  of  250  m,  equivalent  to  the  pre¬ 
vious  two  final  formations.  This  reconfiguration  was  the  largest  performed,  requiring  16 
orbital  periods  for  the  control  time  to  keep  from  exceeding  the  maximum  input  allowed. 
Even  with  the  large  reconfiguration  time,  the  control  input  for  this  reconfiguration  was  still 
relatively  large,  as  shown  in  Figure  4.35,  requiring  a  maximum  ballistic  coefficient  input 
of  0.09  m2 /kg.  The  fuel  budget  resulted  in  a  total  Av  fuel  budget  of  0.10148m/sec  for 
satellite  1  and  0. 13653  mj  sec  for  satellite  2.  The  resulting  reconfiguration,  which  is  plotted 
in  the  radial/in-track  plane  for  both  the  fixed  and  falling  reference  frames  in  Figures  4.36 
and  4.37,  had  some  interesting  qualities.  Satellite  1  took  advantage  of  its  favorable  initial 
conditions,  slowly  increasing  its  in-track  separation  to  allow  satellite  2  time  to  overcome 
its  unfavorable  initial  position.  But  both  satellites  also  reduced  their  radial  oscillations 
as  satellite  1  increased  its  semimajor  axis  and  satellite  2  decreased  its  semimajor  axis, 
allowing  the  separation  distance  between  the  two  satellites  to  decrease.  Then,  as  the  in¬ 
track  distance  continued  to  decrease  toward  the  required  separation,  both  satellites  began 
increasing  their  radial  oscillations  again  to  meet  the  final  formation  conditions  as  their 
average  radial  position  moved  back  to  zero.  These  time  history  plots  of  the  radial  and 
in-track  positions  are  shown  in  Figure  4.38. 

Finally,  as  a  check  to  ensure  the  desired  formation  resulted  from  the  control,  the 
resulting  state  was  propagated  over  20  orbital  periods.  A  plot  of  the  radial/in-track  plane 
motion  over  the  20  periods  is  shown  in  Figure  4.39.  The  two  satellites  are  tracing  out  the 
same  elliptical  formation  with  a  semimajor  axis  of  250  m  as  a  result  of  the  control  input, 
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and  the  formation  is  not  drifting  apart.  This  negligible  drift  is  confirmed  in  the  radial  and 
in-track  position  deviation  plots  in  Figure  4.40.  These  plots  confirm  that  the  satellites  are 
in  the  proper  formation,  although  the  resulting  eccentricities  slightly  differ  from  the  ones 
prescribed,  resulting  in  a  5  mm  radial  position  deviation.  Also,  a  sub-millimeter  difference 
between  the  two  semimajor  axes  results,  leading  to  a  negative  drift  between  the  satellites 
that  brings  the  satellite  1  cm  closer  over  the  20  orbital  periods.  Again,  these  plots  confirm 
the  successful  reconfiguration,  meeting  the  prescribed  constraints  and  achieving  desired 
formations. 

Table  4.6  Payload  Separation  to  Elliptical  Reconfiguration  Details,  6800  km 


Reference  Orbital  Radius 

6800  km 

Initial  Conditions 

xi  =  [—10m,  20 m,  — 0.04m/s,  0.005m/s] 
X2  =  [10  m,  —20  m,  0.04  m/s,  —  0.005  m/s] 

Final  Formation 

Elliptical:  250  m  semimajor  axis 

Total  Time  for 

16  Orbital  Periods 

Reconfiguration 

(24  hr,  48  min ,  8  sec ) 

Maneuver 

Aui  =  0.10148  ml  sec 

Cost 

Av2  =  0.13653  ml  sec 

Altitude  Loss 

197.70  m 

Control  Cost 

21.08 

Constraint  Cost 

1 

03 

OO 

Figure  4.35 


Control  Input:  Payload  Separation  to  Elliptical,  6800  km 
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Figure  4.36  Fixed  Frame  Plot:  Payload  Separation  to  Elliptical,  6800  km 


Figure  4.37  Falling  Frame  Plot:  Payload  Separation  to  Elliptical,  6800  km 


Figure  4.38  Reconfiguration  Time  History  Plot:  Payload  Separation  to  Elliptical, 
6800  km 
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Figure  4.39  Post  Transfer  Plot:  Payload  Separation  to  Elliptical,  6800  km 


Figure  4.40  Post  Transfer  Time  History:  Payload  Separation  to  Elliptical,  6800  km 


V.  Conclusions  and  Recommendations 


A  linear  quadratic  terminal  controller  was  used  to  optimally  reconfigure  a  satellite  for¬ 
mation  using  atmospheric  drag  actuated  control  while  minimizing  the  loss  of  energy  of 
each  satellite  in  the  formation.  Since  the  energy  is  directly  related  to  orbital  altitude,  this 
control  law  effectively  minimized  the  loss  of  altitude  as  a  result  of  the  reconfiguration. 
Reconfigurations  to  in-plane  and  elliptical  formations  were  performed  at  an  orbital  alti¬ 
tude  of  6800  km  from  three  separate  initial  formations:  an  in-plane  formation,  an  elliptical 
formation,  and  a  post-payload  separation  state  whose  initial  conditions  were  selected  to  re¬ 
semble  the  state  of  each  spacecraft  shortly  after  separating  from  the  launch  vehicle.  These 
six  reconfigurations  went  from  small  to  large  formations,  but  Appendix  A  presents  results 
for  the  reverse  reconfigurations,  excluding  the  post-payload  separation  initial  state.  Addi¬ 
tional  results  were  compiled  for  reconfigurations  from  in-plane  and  elliptical  formations  at 
an  orbital  altitude  of  7000  km.  These  results  are  presented  in  Appendix  B.  In  each  case, 
the  altitude  loss  and  the  total  control  effort,  measured  in  terms  of  the  conventional  Av 
budget,  were  recorded  as  performance  measures  for  the  reconfiguration.  Upon  completion 
of  the  reconfiguration,  the  final  formation  state  was  propagated  forward  over  20  orbital 
periods  to  ensure  the  desired  formation  and  constraints  were  satisfied. 

The  results  of  this  study  show  that  the  radial  and  in-track  motion  of  a  satellite  for¬ 
mation  can  be  effectively  controlled  using  an  atmospheric  drag  actuated  control  scheme, 
whereas  the  cross-track  motion  is  uncontrollable,  as  shown  in  Chapter  III.  Additionally, 
minimizing  the  loss  of  energy  conserves  the  loss  of  altitude  during  the  reconfiguration. 
The  desired  formations  were  met  for  each  case,  although  the  time  requirement  varied  de¬ 
pending  upon  the  initial  and  final  conditions.  The  change  in  orbital  eccentricity  proved  to 
dictate  the  control  time  for  a  given  altitude,  and  the  time  requirement  increased  substan¬ 
tially  with  altitude.  Further  increasing  the  allotted  time  for  the  reconfiguration  tended 
to  tighten  the  final  constraints  for  each  transfer  performed.  The  control  law  also  proved 
sufficient  in  limiting  the  drift,  even  though  the  terminal  constraints  were  kept  soft  in  the 
optimization  routine.  The  maximum  post-reconfiguration  drift  due  to  the  final  conditions 
of  the  optimization  routine  was  only  8  cm  over  the  20  orbital  periods. 
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Of  all  the  reconfigurations  performed,  the  in-plane  to  in-plane  transfer  proved  the 
most  efficient.  For  the  6800  km  reference  orbit,  the  formation  lost  only  18.41m  of  altitude 
over  the  4  orbital  period  time  of  reconfiguration.  The  strategy  for  this  transfer  provides 
intuition  into  its  simplicity.  Satellite  1,  the  leading  satellite  at  transfer  initiation,  deployed 
its  panels  first,  letting  the  drag  acceleration  decrease  its  altitude  to  the  final  altitude  before 
zeroing  the  panel  orientation.  The  lower  altitude  allowed  the  separation  distance  between 
the  satellites  to  increase  to  the  desired  500  m  value  as  satellite  2  maneuvered  down  to 
the  final  altitude.  For  the  in-plane  to  in-plane  transfer  presented  in  Section  A. 1.1  from 
an  initial  separation  distance  of  500  m  to  a  final  separation  of  50  m,  which  is  the  reverse 
reconfiguration  of  the  transfer  shown  in  Section  4.2.1,  the  control  strategy  was  exactly  the 
same,  just  applied  to  the  opposite  satellite.  In  this  case,  satellite  2  deployed  its  panels 
first  to  close  the  separation  distance  between  the  two  satellites  before  satellite  1  initiated 
control,  resulting  in  the  same  altitude  loss  and  Av  budget.  More  importantly,  the  control 
effort  was  simplified  since  an  initial  radial  oscillation  from  an  orbital  eccentricity  difference 
was  not  present,  causing  the  controller  to  overcome  a  large  eccentricity. 

This  large  eccentricity  dampening  can  be  seen  in  Section  A.  1.2  where  an  initial 
elliptical  formation  for  a  250  m  semimajor  axis  is  reconfigured  to  an  in-plane  formation 
with  a  50  m  separation  distance.  The  initial  orbit  of  each  satellite  has  a  nonzero  eccentricity 
that  must  be  zeroed  as  the  satellites  approach  the  desired  50  m  separation  in  the  same, 
circular  orbit.  A  substantial  control  effort  is  required  to  circularize  the  orbit,  resulting  in 
an  altitude  loss  of  approximately  150  m.  However,  this  substantial  control  effort  is  required 
due  to  the  size  of  the  initial  eccentricity.  In  Section  4.2.2,  the  initial  formation  is  elliptical 
with  only  a  25  m  semimajor  axis,  resulting  in  a  smaller  eccentricity  than  the  previous. 
The  control  effort  to  overcome  this  radial  oscillation  was  not  nearly  as  substantial,  and  the 
controller  was  actually  able  to  take  advantage  of  the  oscillation  to  increase  the  separation 
distance  of  the  satellites  while  circularizing  the  orbits.  This  transfer  resulted  in  a  total 
altitude  loss  of  only  20  m. 

The  higher  altitude  orbits  posed  the  most  significant  problem  for  the  control  effort. 
Due  to  the  substantially  lower  atmospheric  density  when  the  altitude  was  increased  200  m, 
the  controller  required  much  more  time  to  perform  the  reconfiguration.  A  reconfiguration 
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from  an  elliptical  formation  with  a  25  m  semimajor  axis  to  an  in-plane  formation  with  a 
separation  distance  of  250  m  required  21  orbital  periods,  or  roughly  34  hours,  to  complete 
the  transfer.  This  control  time  increased  when  an  eccentricity  had  to  be  imposed  on  a 
circular  orbit.  A  reconfiguration  from  an  in-plane  formation  with  a  separation  of  50  m 
to  an  elliptical  formation  with  a  semimajor  axis  of  75  m  required  54  orbital  periods  to 
perform  the  transfer,  a  control  time  of  over  87  hours.  Although  the  control  time  increased 
substantially  for  these  test  cases,  the  altitude  loss  remained  within  the  bounds  of  the 
previous  reconfigurations  since  this  value  is  dependent  only  on  the  overall  control  input 
applied.  The  time  requirement  could  be  reduced  if  the  maximum  differential  ballistic 
coefficient  achievable  for  the  formation  was  increased.  This  study  used  a  conservative 
maximum  value  of  0.1  in2 /kg  for  the  differential  ballistic  coefficient  to  provide  a  middle 
point  for  drag  actuated  control.  But  this  value  could  be  increased  through  an  efficient 
spacecraft  design,  resulting  in  decreased  reconfiguration  times  at  all  altitudes. 

Future  research  efforts  should  focus  on  further  validating  drag  actuated  control  and 
its  potential  cost  savings.  This  study  was  primarily  a  proof  of  concept  study.  As  such,  a 
simple  model  was  applied  to  prove  that  reasonable  results  are  obtainable  for  the  prescribed 
control  approach.  To  further  validate  both  the  effectiveness  and  usefulness  of  this  control 
technique,  high  fidelity  atmospheric  and  gravity  models  should  be  applied  with  a  control 
technique  to  account  for  these  variations  in  the  motion.  A  study  that  encompasses  a  full 
range  of  inclinations,  and  perhaps  even  accounts  for  different  reference  orbit  eccentricities, 
could  help  expand  this  concept. 

With  a  full  dynamic  model,  formation  degradation  is  likely  to  occur  on  some  time 
scale.  This  study  argued  that  through  the  proper  formation  design  and  selection  of  or¬ 
bital  elements,  this  degradation  could  be  effectively  minimized.  A  formation  keeping  cost 
comparison  between  drag  actuated  controlled  formations  and  conventionally  controlled 
formations  using  high  impulse  thrusters  could  further  expand  this  concept.  Past  research 
shows  that  a  low-thrust,  continuous  input  control  law  will  likely  result  in  a  higher  Av  bud¬ 
get  than  high  impulse,  finite  control  inputs,  but  these  costs  cannot  be  directly  compared 
since  they  relate  to  entirely  different  costs;  drag  actuated  control  costs  relate  to  altitude 
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while  conventional  control  costs  are  associated  with  onboard  expendables.  So  a  proper 
characterization  of  the  costs  would  involve  the  time  on-orbit  with  available  control. 

Another  interesting  study  would  compare  the  potential  fuel  savings  formation  flying 
with  drag  actuated  control  can  provide  over  formation  flying  using  conventional  thruster. 
When  flying  with  drag,  the  primary  source  of  fuel  is  provided  by  altitude.  So  the  remaining 
lifetime  of  the  satellite  is  characterized  by  its  current  altitude.  For  formation  flying  with 
conventional  thrusters,  the  fuel  onboard  the  satellite  is  usually  depleted  before  the  satellite 
falls  out  of  orbit,  allowing  a  characterization  of  its  lifetime  based  on  the  remaining  fuel 
onboard.  Since  the  drag  actuated  control  satellites  do  not  require  much  fuel  onboard,  they 
would  be  lighter  in  weight  at  launch,  reducing  the  cost  of  orbital  injection.  Moreover,  the 
overall  satellite  could  be  smaller  in  design  since  the  bulky  fuel  tanks  would  not  be  needed, 
allowing  a  single  launch  vehicle  to  lift  more  satellites  into  orbit.  However,  to  control  the 
formation,  the  satellites  are  sacrificing  altitude,  which  may  require  the  launch  vehicle  to 
place  the  satellites  into  a  higher  orbit  to  optimize  the  formation  time  on  orbit.  So  a 
comparison  of  the  costs  required  to  place  a  drag  actuated  controlled  satellite  formation  in 
a  higher  orbit  versus  placing  a  conventionally  control  formation  with  a  greater  total  mass 
would  provide  insight  into  the  overall  cost  savings  of  drag  controlled  system.  Additionally, 
the  time  on  orbit  for  each  scenario  can  be  derived  from  a  common  mission  schedule  requiring 
an  equal  number  of  reconfigurations  and  offering  the  same  formation  keeping  constraints 
as  a  separate  cost  analysis. 

Overall,  the  results  of  the  study  confirmed  that  a  satellite  formation  can  be  effectively 
controlled  in  the  radial/in-track  plane  using  atmospheric  drag  actuated  control  while  the 
overall  energy,  which  is  directly  tied  to  the  altitude  of  the  formation,  is  conserved.  The 
only  drawback  to  such  a  control  effort  is  the  time  requirement  for  the  reconfigurations. 
However,  if  the  reconfiguration  time  is  not  an  issue,  atmospheric  drag  proves  to  be  an  ideal 
and  effective  control  source  for  satellite  formations  in  low-Earth  orbit. 
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Appendix  A.  Additional  Reconfigurations  at  6800  km 
A.l  Reconfiguration  to  an  In-Plane  Formation 

This  section  contains  the  results  for  reconfigurations  into  a  final  in-plane  formation 
with  a  separation  distance  of  500  m  from  two  separate  initial  formations:  an  in-plane 
formation  with  a  separation  distance  of  50  m  and  an  elliptical  formation  with  a  semimajor 
axis  of  25  m  such  that  the  satellites  are  50  m  apart  in  the  in-track  position.  The  satellite 
positions  and  motion  are  with  respect  to  a  6800  km  reference  orbit. 

A.  1.1  In- Plane  to  In- Plane.  The  following  figures  show  the  results  of  a  recon¬ 
figuration  from  an  in-plane  formation  with  a  separation  distance  of  500  m  to  a  smaller 
in-plane  formation  with  a  separation  distance  of  50  m: 

Table  A.l  Large  In-Plane  to  Small  In-Plane  Reconfiguration  Details,  6800  km 


Reference  Orbital  Radius 

6800  km 

Initial  Formation 

In-Plane:  500  m  separation 

Final  Formation 

In-Plane:  50  m  separation 

Total  Time  for 

4  Orbital  Periods 

Reconfiguration 

(6  hr,  12  min ,  2  sec ) 

Maneuver 

Avi  =  0.01048  ml  sec 

Cost 

Av2  =  0.01048  ml  sec 

Altitude  Loss 

18.41  m 

Control  Cost 

0.65 

Constraint  Cost 

le— 5 

Figure  A.l  Control  Input:  Large  In-Plane  to  Small  In-Plane,  6800  km 
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Figure  A. 2  Fixed  Frame  Plot:  Large  In-Plane  to  Small  In-Plane,  6800  km 


Figure  A. 3  Falling  Frame  Plot:  Large  In-Plane  to  Small  In-Plane,  6800  km 


Figure  A. 4  Reconfiguration  Time  History  Plot:  Large  In-Plane  to  Small  In-Plane, 
6800  km 
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Figure  A. 5  Post  Transfer  Plot:  Large  In-Plane  to  Small  In-Plane,  6800 /cm 


Figure  A. 6  Post  Transfer  Time  History:  Large  In-Plane  to  Small  In-Plane,  6800  km 


A.  1.2  Elliptical  to  In- Plane.  The  following  figures  show  the  results  of  a  reconfigu¬ 
ration  from  an  elliptical  formation  with  a  semimajor  axis  of  250  m  to  an  in-plane  formation 
with  a  separation  distance  of  50  m.  Note  that  the  post-transfer  plot  shows  the  satellites 
drifting  over  time.  However,  as  shown  in  the  in-track  drift  plot,  the  satellites  are  actually 
drifting  in  the  same  direction  such  that  their  total  separation  after  20  orbital  periods  is 
only  2.2  cm. 

Table  A. 2  Large  Elliptical  to  Small  In-Plane  Reconfiguration  Details,  6800  km 


Reference  Orbital  Radius 

6800  km 

Initial  Formation 

Elliptical:  250  m  semimajor  axis 

Final  Formation 

In-Plane:  50  m  separation 

Total  Time  for 

10  Orbital  Periods 

Reconfiguration 

(15  hr,  30  min,  5  sec) 

Maneuver 

A  ci  =  0.08964  ml  sec 

Cost 

Av2  =  0.08964  ml  sec 

Altitude  Loss 

148.92  m 

Control  Cost 

18.32 

Constraint  Cost 

le-5 

Time  (Orbital  Periods) 


Figure  A. 7  Control  Input:  Large  Elliptical  to  Small  In-Plane,  6800  km 


A-4 


Figure  A. 8  Fixed  Frame  Plot:  Large  Elliptical  to  Small  In-Plane,  6800  km 


Figure  A. 9  Falling  Frame  Plot:  Large  Elliptical  to  Small  In-Plane,  6800  km 


Figure  A.  10  Reconfiguration  Time  History  Plot:  Large  Elliptical  to  Small  In-Plane, 
6800  km 
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Figure  A. 11  Post  Transfer  Plot:  Large  Elliptical  to  Small  In-Plane,  6800  km 


Figure  A. 12  Post  Transfer  Time  History:  Large  Elliptical  to  Small  In-Plane,  6800  km 


A. 2  Reconfiguration  to  an  Elliptical  Formation 

This  section  contains  the  results  for  reconfigurations  into  a  final  elliptical  formation 
with  a  semimajor  axis  of  25  m,  such  that  the  satellites  are  50  m  apart  when  aligned  in 
the  in-track  direction,  from  two  separate  initial  formations:  an  in-plane  formation  with  a 
separation  distance  of  500  m  and  an  elliptical  formation  with  a  semimajor  axis  of  250  m. 
The  satellite  positions  and  motion  are  with  respect  to  a  6800  km  reference  orbit. 

A. 2.1  In-Plane  to  Elliptical.  The  following  figures  show  the  results  of  a  recon¬ 
figuration  from  an  in-plane  formation  with  a  separation  distance  of  500  m.  to  an  elliptical 
formation  with  a  semimajor  axis  of  25  m: 

Table  A. 3  Large  In-Plane  to  Small  Elliptical  Reconfiguration  Details,  6800  km 


Reference  Orbital  Radius 

6800  km 

Initial  Formation 

In-Plane:  500  m  separation 

Final  Formation 

Elliptical:  25  m  semimajor  axis 

Total  Time  for 

6  Orbital  Periods 

Reconfiguration 

(9  hr ,  18  min ,  3  sec) 

Maneuver 

Am  =  0.01151  m/sec 

Cost 

Av2  =  0.01151  m/sec 

Altitude  Loss 

19.12  m 

Control  Cost 

0.61 

Constraint  Cost 

4e— 6 

Figure  A. 13  Control  Input:  Large  In-Plane  to  Small  Elliptical,  6800  km 
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Figure  A. 14  Fixed  Frame  Plot:  Large  In-Plane  to  Small  Elliptical,  6800  km 


Figure  A. 15  Falling  Frame  Plot:  Large  In-Plane  to  Small  Elliptical,  6800  km 


Figure  A.  16  Reconfiguration  Time  History  Plot:  Large  In-Plane  to  Small  Elliptical, 

6800  km 
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Figure  A. 17  Post  Transfer  Plot:  Large  In-Plane  to  Small  Elliptical,  6800  km 


Figure  A.  18  Post  Transfer  Time  History:  Large  In-Plane  to  Small  Elliptical,  6800  km 


A. 2.2  Elliptical  to  Elliptical.  The  following  figures  show  the  results  of  a  recon¬ 
figuration  from  an  elliptical  formation  with  a  250  m  semimajor  axis  to  a  smaller  elliptical 
formation  with  a  25  m  semimajor  axis.  Note  that  the  post  transfer  plot  shows  a  drift  in 
the  formation  over  time.  However,  as  shown  in  the  in-track  deviation  plot  in  Figure  A. 24, 
the  two  satellites  are  drifting  together,  resulting  in  a  total  positive  drift  of  only  3  cm. 

Table  A. 4  Large  Elliptical  to  Small  Elliptical  Reconfiguration  Details,  6800  km . 


Reference  Orbital  Radius 

6800  km 

Initial  Formation 

Elliptical:  250  m  semimajor  axis 

Final  Formation 

Elliptical:  25  m  semimajor  axis 

Total  Time  for 

8  Orbital  Periods 

Reconfiguration 

(12  hr,  24  'min ,  4  sec) 

Maneuver 

Avi  =  0.08069  ml  sec 

Cost 

Av2  =  0.08069  ml  see 

Altitude  Loss 

134.05  m 

Control  Cost 

18.60 

Constraint  Cost 

2e— 5 

Time  (Orbital  Periods) 


Figure  A.  19  Control  Input:  Large  Elliptical  to  Small  Elliptical,  6800  km 
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Figure  A. 20  Fixed  Frame  Plot:  Large  Elliptical  to  Small  Elliptical,  6800  km 
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Figure  A. 21  Falling  Frame  Plot:  Large  Elliptical  to  Small  Elliptical,  6800  km 


Figure  A. 22  Reconfiguration  Time  History  Plot:  Large  Elliptical  to  Small  Elliptical, 

6800  /cm 


Figure  A. 23  Post  Transfer  Plot:  Large  Elliptical  to  Small  Elliptical,  6800  km 
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Figure  A. 24  Post  Transfer  Time  History:  Large  Elliptical  to  Small  Elliptical,  6800  km 


Appendix  B.  Reconfigurations  at  7000  km 
B.l  Reconfiguration  to  an  In- Plane  Formation 

This  section  contains  the  results  for  reconfigurations  into  a  final  in-plane  formation 
with  a  separation  distance  of  250  m  from  two  separate  initial  formations:  an  in-plane 
formation  with  a  separation  distance  of  50  m  and  an  elliptical  formation  with  a  semimajor 
axis  of  25  m  such  that  the  satellites  are  50  m  apart  in  the  in-track  position.  The  satellite 
positions  and  motion  are  with  respect  to  a  7000  km  reference  orbit. 

B.1.1  In- Plane  to  In-Plane.  The  following  figures  show  the  results  of  a  reconfig¬ 
uration  from  an  in-plane  formation  with  a  separation  distance  of  50  m  to  another  in-plane 
formation  with  a  separation  distance  of  250  m  at  a  reference  orbital  radius  of  7000  km: 

Table  B.l  In-Plane  to  In-Plane  Reconfiguration  Details,  7000 km. 


Reference  Orbital  Radius 

7000  km 

Initial  Formation 

In-Plane:  50  m  separation 

Final  Formation 

In-Plane:  250  m  separation 

Total  Time  for 

7  Orbital  Periods 

Reconfiguration 

(11  hr,  19  min,  59  sec) 

Maneuver 

Avi  =  0.00246  ml  sec 

Cost 

Av2  =  0.00246  ml  sec 

Altitude  Loss 

4.55  m 

Control  Cost 

0.02 

Constraint  Cost 

2e— 7 

Figure  B.l  Control  Input:  In-Plane  to  In-Plane,  7000 km 
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Figure  B.2  Fixed  Frame  Plot:  In-Plane  to  In-Plane,  7000  km 
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Figure  B.3  Falling  Frame  Plot:  In-Plane  to  In-Plane,  7000 /cm 


Figure  B.4  Reconfiguration  Time  History  Plot:  In-Plane  to  In-Plane,  7000  km 


B.1.2  Elliptical  to  In-Plane.  The  following  figures  show  the  results  of  a  reconfig¬ 
uration  from  an  elliptical  formation  with  a  25  m  semimajor  axis  to  an  in-plane  formation 
with  a  separation  distance  of  250  m  at  a  reference  orbital  radius  of  7000  km 

Table  B.2  Elliptical  to  In-Plane  Reconfiguration  Details,  7000  km 


Reference  Orbital  Radius 

7000  km 

Initial  Formation 

Elliptical:  25  m  semimajor  axis 

Final  Formation 

In-Plane:  250  m  separation 

Total  Time  for 

21  Orbital  Periods 

Reconfiguration 

(33  hr,  59  min,  59  sec ) 

Maneuver 

Aui  =  0.00863  ml  sec 

Cost 

Av2  =  0.00863  ml  sec 

Altitude  Loss 

15.92  m 

Control  Cost 

0.08 

Constraint  Cost 

le-8 

Figure  B.7  Control  Input:  Elliptical  to  In-Plane,  7000  km 


Figure  B.8  Fixed  Frame  Plot:  Elliptical  to  In-Plane,  7000  km 
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B.2  Reconfiguration  to  an  Elliptical  Formation 

This  section  contains  the  results  for  reconfigurations  into  a  final  elliptical  formation 
with  a  semimajor  axis  of  75  m,  such  that  the  satellites  are  150  m  apart  when  aligned  in 
the  in-track  direction,  from  two  separate  initial  formations:  an  in-plane  formation  with  a 
separation  distance  of  50  nn  and  an  elliptical  formation  with  a  semimajor  axis  of  25  m.  The 
satellite  positions  and  motion  are  with  respect  to  a  7000  km.  reference  orbit. 

B.2.1  In-Plane  to  Elliptical.  The  following  figures  show  the  results  of  a  recon¬ 
figuration  from  an  in-plane  formation  with  a  separation  distance  of  50  m  to  an  elliptical 
formation  with  a  75  m  semimajor  axis  at  a  reference  orbital  radius  of  7000  km: 

Table  B.3  In-Plane  to  Elliptical  Reconfiguration  Details,  7000  km 


Reference  Orbital  Radius 

7000  km 

Initial  Formation 

In-Plane:  50  m  separation 

Final  Formation 

Elliptical:  75  m  semimajor  axis 

Total  Time  for 

54  Orbital  Periods 

Reconfiguration 

(87  hr ,  25  min ,  40  sec ) 

Maneuver 

Am  =  0.02574  m/sec 

Cost 

Av2  =  0.02574  m/sec 

Altitude  Loss 

44.52  m 

Control  Cost 

0.26 

Constraint  Cost 

le— 8 

Time  (Orbital  Periods) 


Figure  B.13  Control  Input:  In-Plane  to  Elliptical,  7000  km 


B-7 


B.2.2  Elliptical  to  Elliptical.  The  following  figures  show  the  results  of  a  recon¬ 
figuration  from  an  elliptical  formation  with  a  25  m  semimajor  axis  to  another  elliptical 
formation  with  a  75  m  semimajor  axis  at  a  reference  orbital  radius  of  7000  km: 


Table  B.4  Elliptical  to  Elliptical  Reconfiguration  Details 


Reference  Orbital  Radius 

7000  km 

Initial  Formation 

Elliptical:  25  m  semimajor  axis 

Final  Formation 

Elliptical:  75  m  semimajor  axis 

Total  Time  for 

36  Orbital  Periods 

Reconfiguration 

(58  hr,  17  min,  7  sec ) 

Maneuver 

Av\  =  0.01716  ml  sec 

Cost 

Av2  =  0.01716  ml  sec 

Altitude  Loss 

29.68  m 

Control  Cost 

0.17 

Constraint  Cost 

le-8 

Figure  B.19  Control  Input:  Elliptical  to  Elliptical,  7000  km 


Figure  B.20  Fixed  Frame  Plot:  Elliptical  to  Elliptical,  7000  km 
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Figure  B.21  Falling  Frame  Plot:  Elliptical  to  Elliptical,  7000  km 


Figure  B.22  Reconfiguration  Time  History  Plot:  Elliptical  to  Elliptical,  7000  km 
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Appendix  C.  Source  Code 

Shown  below  is  the  Matlab®  source  code  used  to  generate  the  results  presented  within 
this  study. 

•/. -  TERMINAL_CONTROLLER  - 

'/,  Terminal  Controller:  Main  piece  of  code  used  to  determine  the  optimal 
'/,  reconfiguration  of  the  user-specified,  two-satellite  formations. 

'/,  The  user  specifies  whether  to  drive  formation  to  an  in-plane,  elliptical 
'/,  (2x1)  ,  or  a  formation  with  zero  in-track  drift  but  no  constraints  on  the 
'/,  radial  oscillations,  from  some  given  initial  formation  or  some  initial 
'/,  position  for  each  satellite.  Takes  the  final  state  from  the  optimization 
'/,  routine  and  propagates  it  forward  in  time  using  the  mean  motion  for  the 
'/,  new  altitude  to  show  that  the  assumption  of  a  time-invariant  system 
•/.  holds. 

'/,  Code  calls  the  following  subfunctions: 

'/,  -  atmos_exp.m:  exponential  atmospheric  density  model 

'/,  -  cw2_eom.m:  differential  equations  for  the  two-satellite  formation 

'/,  -  ham_eom.m:  Hamiltonian  system  of  equations  for  the  two-satellite 

'/,  formation 

'/,  Written  by:  Blake  Hajovsky 

'/,  —  Original  version:  10  Oct  06 

'/,  — >  Last  update :  15  Feb  07 

function  terminal_controller 
clear  all 

'/,  Define  the  global  variables: 
global  nO  rho  u  vO 
tic  "/,  Start  the  clock 

'/. -  USER  INPUTS  - 

'/,  Define  the  initial  formation: 

'/,  init_form  =  1  =>  in-plane  formation  (no  in-track  drift/radial  osc.) 
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'/,  init_form  =  2  =>  elliptical  (2x1  ellipse;  no  in-track  drift) 

'/,  init_form  =  3  =>  post  payload  separation  or  any  other  user  defined 

'/,  initial  conditions 

init_form  =  2; 

'/,  If  the  in-plane  or  elliptical  formation  was  chosen,  define  the  in-track 
'/,  separation  distance  for  the  initial  formation: 
d_sepO  =  50; 

'/,  If  the  post-payload  separation  formation  or  any  other  formation  was 
'/,  selected  for  the  initial  orbit,  define  the  initial  relative  position  of 
'/,  the  satellites  (meters)  : 

'/,  xO  =  [x  y  z  |  dx  dy  dz]  ’ 

xOl  =  [-10;  20;  0;  -0.04;  0.005;  0];  x02  =  [10;  -20;  0;  0.04; 

-0.005;  0] ; 

'/,  Define  the  final  formation: 

'/,  fig  =  1  =>  in-plane  (no  in-track  drift;  no  radial  oscillation) 

'/,  fig  =  2  =>  elliptical  (2x1  ellipse;  no  in-track  drift) 

'/,  fig  =  3  =>  formation  without  radial  constraint  (w/  no  in-track  drift) 

fig  =  2; 

'/,  Define  the  final  in-track  separation  distance: 
d_sepf  =  500; 

'/,  Define  the  optimization  parameters  and  other  values: 

'/,  Define  the  circular  reference  orbit  radius  (km)  : 
rO  =  6800;  /  km 

'/,  Define  the  length  of  time  over  which  to  optimize  the  system  (orbital 
'/,  periods)  : 
t_final  =  8; 

'/,  Set  the  number  of  steps  for  the  linear  propagation  of  the  controlled 
'/,  system: 

Ns  =  5000; 

'/,  Define  the  length  of  time  to  propagate  the  E0M  after  the  optimization 
'/,  (orbital  periods)  : 
t_prop  =  20; 

'/,  Set  the  weight  on  the  final  constraints  (wed  =  weight  on  drift, 

'/,  wes  =  weight  on  separation  distance)  : 
wed  =  8e8;  wes  =  8e8; 
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'/,  Specify  whether  or  not  to  plot  the  ’fixed’  plane,  holding  the  relative 
'/,  reference  frame  constant  (positive  control  inputs  only)  : 

'/,  NOTE:  Keep  OFF  if  not  being  used  to  speed  up  the  simulation. 

'/,  fixed_frame  =  0  =>  OFF 

'/,  fixed_frame  =  1  =>  ON 

fixed_frame  =  0; 

'/,  Specify  where  to  create  the  folder  that  the  *.eps  files  will  be  sent: 
foldername  =  [’el2el_l’];  mkdir (foldername) 

'/. -  END  USER  INPUTS  - 

'/,  We  should  print  to  the  screen  what  we  have  selected  with  the  input  from 
'/,  above  as  a  check  on  what  we  are  actually  doing.  Do  this  here: 

'/,  First,  print  the  initial  formation: 

if  init_form  ==  1  "/,  an  in-plane  formation  has  been  selected 

disp(sprintf  (’  Initial  Formation:  In-plane  w/  separation  of  "/,  4.2f  m.’,d_sep0)) 
elseif  init_form  ==  2  "/,  an  elliptical  formation  has  been  selected 

disp(sprintf  (’  Initial  Formation:  Elliptical  w/  ’’a’’  =  °/04.2f  m.  ’  ,0.5*d_sep0)) 
if  startl  ==  0,  disp( ’Transfer  begins  from  the  RADIAL  position.’) 
elseif  startl  ==  1,  disp( ’Transfer  begins  from  the  IN-TRACK  position.’) 
else 

'/,  An  incorrect  selection  has  been  made,  so  break  the  code  here  so 
'/,  the  problem  can  be  corrected: 

error(’An  incorrect  selection  has  been  made  for  ’’startl,’’  where  to 
begin  the  transfer.’) 

end 

elseif  init_form  ==  3  "/,  the  user  chose  the  initial  conditions 

disp(sprintf (’ Initial  Formation:  User  defined.’)); 

else 

"/,  An  incorrect  selection  has  been  made,  so  break  the  code  here  so  the 
"/,  problem  can  be  corrected: 

error(’An  incorrect  selection  has  been  made  for  ’ ’ init_f orm, ’ ’  the  initial 
formation. ’ ) 

end 

'/,  Now,  print  the  final  formation: 
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if  fig  ==  1  7.  an  in-plane  formation 

disp(sprintf ( ’Final  Formation:  In-plane  w/  separation  of  7,4.2f  m.  ’  ,d_sepf )  )  ; 
elseif  fig  ==  2  7  an  elliptical  formation 

disp(sprintf  ( ’Final  Formation:  Elliptical  w/  ’’a’’  =  '/,4.2f  m. ’ , 0 . 5*d_sepf ) ) ; 
elseif  fig  ==  3  "/,  formation  w/out  a  radial  constraint 

disp(sprintf ( ’Final  Formation:  In-track  separation  distance  only  constrained 
w/  d  =  "/, 4.2f  m’,d_sepf)); 

else 

7,  An  incorrect  selection  has  been  made,  so  break  the  code  here  so  the 
"/,  problem  can  be  corrected: 

error(’An  incorrect  selection  has  been  made  for  ’’fig,’’  the  final  formation.’) 

end 

'/,  Output  the  reference  orbit  radius  and  controlled  time: 
disp(sprintf  ( ’Reference  Orbital  Radius:  "/,4.0f  km.’,rO)); 
disp(sprintf (’Controlled  over  % 2. Of  orbital  periods ,t_final) ) ; 

7. - BEGIN  MAIN  SCRIPT - 

'/,  Define  the  Earth’s  mean  equatorial  radius  and  the  gravitational 
'/,  parameter: 

Re  =  6378.137;  7.  km 

mu  =  398600.4418;  7.  km~3/sec~2 

7.  Determine  the  mean  motion  of  the  reference  orbit: 
nO  =  sqrt (mu/r0~3) ; 

7.  Determine  the  reference  velocity  (km/sec)  and  convert  to  m/s: 
v  =  sqrt(mu/r0);  v_ref  =  v*1000; 

7.  Assuming  the  atmosphere  rotates  with  the  Earth,  calculate  the  atmospheric 
7.  velocity  in  m/sec  (angular  velocity  calculated  based  on  sidereal  day, 

7.  23h  56m  4.0905s  =  86164.0905  sec): 

w_earth  =  2*pi/ (86164 . 0905) ;  v_atm  =  w_earth*r0*1000 ; 

7.  Now,  determine  the  reference  velocity  with  respect  to  the  atmosphere: 
vO  =  v_ref  -  v_atm; 
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7  Determine  the  orbital  period  of  the  reference  orbit: 
t_period  =  2*pi*sqrt (r0~3/mu) ; 

7  Call  ’atmos_exp.m’  to  determine  the  atmospheric  density  at  the  reference 
7  altitude  (kg/m~3)  : 
rho  =  atmos_exp(rO) ; 

'/. -  INITIAL  FORMATION  SET-UP  - 

7  Determine  the  initial  position  of  the  satellites  based  on  the  initial 
7  formation  selected  and  the  desired  transfer  starting  position: 
if  init_form  ==  1 

"/,  An  in-plane  starting  formation  has  been  selected.  This  particular 
7  formation  has  only  an  in-track  separation  distance  (all  velocities 
"/,  are  zero  and  the  radial/cross-track  components  are  zero)  . 
xOl  =  [0;  0 . 5*d_sep0;  0;  0;  0;  0]; 
x02  =  [0;  -0 . 5*d_sep0 ;  0;  0;  0;  0]; 
elseif  init_form  ==  2 

7  An  elliptical  formation  has  been  selected.  The  initial  conditions 
7  are  dependent  upon  the  point  to  the  start  the  transfer. 

7  This  selection  begins  the  transfer  from  the  in-track  separation 
7  point  (zero  radial).  So,  the  initial  conditions  will  have  an 
7  in-track  separation  and  a  radial  velocity.  Determine  these 
7  values : 

in_sepl  =  0.5*d_sep0; 
in_sep2  =  -0.5*d_sep0; 

7  Now  specify  the  radial  velocities: 
vl  =  0 . 25*d_sep0*n0 ; 
v2  =  -0.25*d_sep0*n0; 

7  Define  the  satellites’  states: 
xOl  =  [0;  in_sepl ;  0;  vl;  0;  0]; 
x02  =  [0;  in_sep2;  0;  v2;  0;  0]; 

elseif  init_form  ==  3,  7  Do  nothing;  initial  conditions  already  specified. 

else,  7  Do  nothing; 

end 
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'/,  Build  the  overall  initial  state  vector: 
xO  =  [xOl;  x02] ; 

'/. - OPTIMIZATION  ROUTINE - 

'/,  Determine  the  desired  separation  rate  to  meet  the  2x1  ellipse  constraint: 
dr  =  0 . 5*d_sepf *n0 ; 

'/,  Define  some  of  the  constants  that  will  be  used  ahead: 
n02  =  n0~2;  v02  =  v0~2;  v03  =  v0~3;  v06  =  v0~6;  rho2  =  rho~2; 

'/,  Define  the  A1  matrix  (6x6)  for  the  state  of  1  satellite: 

A1  =  [0,  0,  0,  1,  0,  0;  0,  0,  0,  0,  1,  0;  0,  0,  0,  0,  0,  1;... 

3*n02,  0,  0,  0,  2*n0,  0;  0,  0,  0,  -2*n0,  0,  0;  0,  0,  -n02,  0,  0,  0]; 

'/,  Now  form  the  overall  A  matrix.  Since  the  mean  motion  corresponds  to  the 
'/,  reference  orbit,  A1  =  A2.  So  form  12x12  A  matrix  using  same  entries: 

A  =  [Al,  zeros(6,6);  zeros(6,6),  Al] ; 

'/,  Define  the  control  matrix  for  each  satellite  (6x1)  : 

B1  =  [0;  0;  0;  0;  -0 . 5*rho*v02 ;  0]; 

'/,  Now  form  the  overall  control  matrix  (12x2): 

B  =  [Bl,  zeros(6,l);  zeros(6,l),  Bl] ; 

'/,  Define  the  control  penalty  matrix  R: 

R  =  [0 . 25*rho2*v06 ,  0;  0,  0 . 25*rho2*v06] ; 

'/,  Define  the  final  time: 
tf  =  t_f inal*t_period; 

'/,  Define  Qf ,  the  final  error  weighting  matrix: 

Qf  =  eye(12)*wcd;  Qf(2,2)  =  wcs;  Qf(8,8)  =  wcs; 

'/,  Define  the  terminal  constraint  matrix  (Mf)  and  the  final  constraint 
'/,  vector  (psi)  based  on  the  formation-type  selected: 
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'/.  fig  =  1  =>  in-plane  (no  in-track  drift;  no  radial  oscillation) 

%  fig  =  2  =>  elliptical  (2x1  ellipse;  no  in-track  drift) 

'/.  fig  =  3  =>  formation  without  radial  constraint  (no  in-track  drift) 

d  =  d_sepf ;  if  fig  ==  1 

7.  Define  Mf ,  the  terminal  constraint  matrix: 

Mf  =  zeros (12 , 12) ;  Mf(l,l)  =  2*n0;  Mf(7,7)  =  -2*n0;  Mf(2,2)  =  1; 

Mf  (8,8)  =  -1;  Mf (4,4)  =  1;  Mf(10,10)  =  -1;  Mf(5,5)  =  1;  Mf(ll,ll)  =  -1 
7.  Define  psi,  the  final  constraint  vector: 
psi  =  zeros(12,l);  psi(2)  =  0.5*d;  psi(8)  =  0.5*d; 
elseif  fig  ==  2 

/  Define  Mf ,  the  terminal  constraint  matrix: 

Mf  =  zeros(12, 12) ;  Mf(l,l)  =  2*n0;  Mf(7,7)  =  -2*n0;  Mf(2,2)  =  1; 

Mf  (8,8)  =  -1;  Mf (4,4)  =  1;  Mf(10,10)  =  -1;  Mf(5,5)  =  1;  Mf(ll,ll)  =  -1 
7.  Define  psi,  the  final  constraint  vector: 

psi  =  zeros(12,l);  psi(2)  =  0.5*d;  psi(8)  =  0.5*d;  psi(4)  =  0.5*dr; 
psi(10)  =  0 . 5*dr ; 
elseif  fig  ==  3 

7.  Define  Mf ,  the  terminal  constraint  matrix: 

Mf  =  zeros(12, 12) ;  Mf(l,l)  =  2*n0;  Mf(7,7)  =  -2*n0;  Mf(2,2)  =  -1; 

Mf  (8,8)  =  1;  Mf (5,5)  =  1;  Mf(ll.ll)  =  -1; 

7.  Define  psi,  the  final  constraint  vector: 
psi  =  zeros(12,l);  psi(2)  =  -0.5*d;  psi(8)  =  -0.5*d; 
else,  error (’ Incorrect  ’’fig’’  designator  for  formation-type!’) 
end 

7.  Define  the  Hamiltonian  Matrix: 

H  =  [A,  -B*inv(R)*B’ ;  zeros (12 , 12)  ,  -A’]; 

7.  Find  the  forward  transition  matrix  of  Euler-Lagrange  eqns,  t=tf  to  t=tO: 
Ph  =  expm(H*tf); 

7.  Partition  the  transition  matrix  to  separate  the  states  and  costates 
7.  components : 

state_length  =  length (xO) ;  nl  =  [1 : state_length] ;  n2  = 

[state_length+l : 2*state_length] ;  tranll  =  Ph(nl,nl);  tranl2  = 

Ph(nl,n2);  tran21  =  Ph(n2,nl);  tran22  =  Ph(n2,n2); 
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'/,  Solve  for  the  initial  costates: 

Sf  =  Mf’*Qf*Mf;  a  =  Mf’*Qf*psi;  lambdaO  = 

inv( (tran22-Sf *tranl2) ) * ( (Sf *tranll-tran21) *xO-a) ; 

'/, -  Linear  Propagation  - 

'/,  Now,  use  the  state  transition  matrix  to  linearly  propagate  the  initial 
'/,  states  (defined)  and  costates  (just  found)  to  the  final  time: 
yO  =  [xO;  lambdaO];  t  =  tf  *  [0 : 1/Ns  :  1]  ; 

'/,  We’re  going  to  multiply  by  R  inverse,  B  transpose,  and  the  time  step 
'/,  quite  a  few  times  within  the  next  loop.  So,  calculate  these  things  in 
'/,  advance  to  speed  things  along. 

R_inv  =  inv(R) ;  B_tran  =  B’;  t_step  =  1/Ns*tf ; 

'/,  Initialize  the  summers : 

costl  =  0;  eneg_cost  =  0;  control_totla  =  0;  control_tot2a  =  0; 

'/,  Begin  the  forward  propagation: 
for  i=l:Ns+l 

y(:,i)  =  expm(H*t (i) ) *y0 ; 
x_lin( : ,i)  =  y(nl,i) ; 
u_lin(:,i)  =  -R_inv*B_tran*y (n2 , i) ; 
costl  =  costl  +  u_lin( : , i) ’ *R*u_lin( : , i) *t_step; 
eneg_cost  =  eneg_cost  +  (0 . 5*rho*v03* (abs (u_lin(l , i) )  + . .  . 
abs (u_lin(2, i) ) ) *t_step) ; 

control_totla  =  control_totla  +  abs (0 . 5*rho*v02*u_lin(l , i) *t_step) ; 
control_tot2a  =  control_tot2a  +  abs (0 . 5*rho*v02*u_lin(2 , i) *t_step) ; 

end 

'/, -  Numerically  Integrated  Solution  - 

'/,  Ideally,  the  analytical  results  for  the  controlled  portion  above  should 
'/,  match  what  we  get  if  the  Hamiltonian  system  is  integrated  forward  in  time 
'/,  with  the  same  inputs  (same  initial  states  and  costates).  So  let’s 
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'/,  numerically  integrate  the  EOM  to  verify  the  results: 


'/,  Set-up  the  parameters  for  ode45: 

xO_ham  =  [xO;  lambdaO] ;  t_fin  =  t_f inal*t_period;  tspan  =  [0 
t_fin] ;  options  =  odeset ( ’RelTol ’ , le-9 , ’ AbsTol ’ , le-9) ; 

'/,  Call  ode45  to  propagate  the  Hamiltonian  system  forward  in  time: 
[t_num,out]  =  ode45 (Sham_eom, tspan, x0_ham, options) ; 

'/,  Separate  the  states  and  costates  from  the  numerical  integration: 
x_num  =  out(:,nl);  lam_num  =  out(:,n2); 

'/,  Grab  and  store  the  final  state  and  costate  from  the  numerical 
'/,  integration: 

[n,m]  =  size (out);  xf_all  =  out(n,:);  xf_num  = 
xf _all (1 : state_length) ;  lamf_num  = 
xf _all (state_length+l : 2*state_length) ; 

'/,  In  the  upcoming  loop,  we’ll  determine  the  control  input  at  each  time 
'/,  step.  To  do  this,  we  need  the  R  and  B  matrices.  So  form  these  outside 
'/,  the  loop  and  then  call  from  within: 

'/,  Recall:  rho2  =  rho~2 

•/.  v06  =  v0“6 

R1  =  [0 . 25*rho2*v06] ;  R2  =  R1 ; 

inv_Rl  =  inv(Rl);  "/,  inverse  of  R1 

inv_R2  =  inv(R2)  ;  "/,  inverse  of  R2 

'/,  Define  the  control  matrix  for  each  satellite  (6x1)  : 

'/,  Recall:  v02  =  v0~2 

B1  =  [0;  0;  0;  0;  -0 . 5*rho* (v02) ;  0];  B2  =  Bl; 

Bl_tran  =  Bl’;  7,  transpose  of  Bl 

B2_tran  =  B2’;  7,  transpose  of  B2 

7,  Determine  the  control  input  at  each  time  step: 

for  jj  =  l:n 

%  Grab  the  costates : 
rowl  =  out (j j , : ) ; 
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coni  =  rowl(13:18); 
con2  =  rowl(19:24); 
lam3(j  j  ,  :  )  =  coni ; 
lam4(jj,:)  =  con2; 

7.  Using  R,  B,  and  the  costate,  determine  the  control  input  to  each  sat: 
u_numl(jj)  =  -inv_Rl*Bl_tran*conl ’ ; 
u_num2(jj)  =  -inv_R2*B2_tran*con2’ ; 

end 

7,  Determine  how  close  we  are  to  meeting  the  terminal  constraints  for  both 
7.  the  linear  propagation  and  the  numerical  integration  (output  results 
7.  later)  : 

termc_drift_lin  =  xf _lin(5) -xf _lin(ll)+2*n0* (xf _lin(l) -xf _lin(7) ) ; 
termc_distance_lin  =  xf _lin(8) -xf _lin(2) ;  termc_drif t_num  = 
xf _num(5) -xf _num(ll)+2*n0* (xf _num(l) -xf _num(7) ) ; 
termc_distance_num  =  xf _num(8) -xf _num(2) ; 

7, -  POST  OPTIMIZATION  ANALYSIS  - 

7.  We  now  have  an  optimal  solution  propagated  from  tO  to  tf .  Let’s  see  how 
7.  well  it  actually  performed. 

7.  Given  this  optimal  input,  we  now  want  to  plot  the  controlled  orbit  in  the 
7.  fixed  reference  frame.  Set-up  the  parameters  here: 
xO_new  =  xO;  tspan  =  [0  t_step] ;  t_count  =  0;  eneg_cost2  =  0; 
control_totl  =  0;  control_tot2  =  0;  for  jj  =  1 : length (u_lin) 
for  i  =  1:2 

if  u_lin(i,jj)  <=  0,  u_lin(i,jj)  =  0; 
else,  u_lin(i,jj)  =  2*u_lin(i , j j ) ; 
end 

end 

7.  This  calculates  the  altitude  loss  based  on  the  new  control  input: 
eneg_cost2  =  eneg_cost2  +  0 . 5*rho*v03*u_lin(l , j j ) *t_step ; 

7.  Now  calculate  the  control  cost  on  each  sat  based  on  the  new  control 
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/  input : 

control_totl  =  control_totl  +  abs (0 . 5*rho*v02*u_lin(l , j j ) *t_step) ; 
control_tot2  =  control_tot2  +  abs (0 . 5*rho*v02*u_lin(2 , j j)*t_step) ; 

7.  If  we  decided  to  check  ourselves  with  the  new  input ,  we  should  run 
7.  this  next  part  (specify  in  the  ’User  Inputs’  at  the  beginning  of  the 
%  script).  The  note  below  the  ’if’  statement  has  a  good  explanation  of 
%  this  loop’s  purpose: 
if  fixed_frame  ==  1 

7,  With  this  ’new’  control  input,  we  want  to  propagate  the  EOM 
7.  forward  again  to  ensure  we  get  the  same  result: 
u  =  [u_lin(l , j j ) ;  u_lin(2 , j j )] ; 

[t_n,x_n]  =  ode45(@cw2_eom, tspan,xO_new, options) ; 

7.  Now,  grab  the  final  values  from  the  above  integration: 
nm  =  length (x_n) ; 
x0_new  =  x_n(nm,:); 
x_new(jj,:)  =  x0_new; 

7.  Increment  the  time  count  and  set  the  new  tspan: 
t_count  =  t_count+t_step ; 
t_span  =  [t_count-t_step,  t_count] ; 
else,  %  Do  nothing 
end 

end 

7.  Calculate  the  cost  of  the  terminal  constraints : 
err_diff  =  Mf *x_lin( : ,length(x_lin) ) -psi ;  cost2  = 
err_dif f ’ *Qf *err_dif f ; 

7.  Calculate  the  change  in  semi-major  axis  as  a  result  of  the  maneuver: 

EnergyO  =  - (mu*1000~3) /(2*r0*1000)  ;  7,  meters~2  /  sec~2 

Energyf2  =  EnergyO  -  0 . 5*eneg_cost ;  7.  meters~2  /  sec~2 

rf2  =  -  (mu*1000~3)  /  (2*Energyf  2)  ;  7,  meters 

altitude_loss2  =  r0*1000  -  rf2; 

7.  Now  propagate  the  final  state  from  the  optimization  forward  over  a  user 
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'/,  defined  time  period  to  see  the  orbit  post-optimization  with  zero  control. 
'/,  Set  up  the  parameters  for  ode45.m: 

xO_new  =  xf _num(l : 12) ;  t_fin2  =  t_prop*t_period;  tspan  =  [0 
t_fin2] ;  options  =  odeset ( ’RelTol ’ , le-9 , ’ AbsTol ’ , le-9) ; 

'/,  Calculate  a  new  mean  motion  from  the  new  altitude  to  use  for  the 
'/,  post-optimization  propagation: 

rf_km  =  rf2*10~-3;  nf  =  sqrt (mu/rf _km~3) ;  delta_n  =  nf-nO; 
percent_change_n  =  delta_n/n0*100 ;  nO  =  nf ; 

'/,  Calculate  a  new  satellite  vel  wrt  the  rotating  atmosphere  from  the  new 
'/,  altitude: 

vf  =  sqrt (mu/rf _km) ;  "/,  km/sec 

vf_ref  =  vf*1000;  "/,  m/sec 

vf_atm  =  w_earth*rf _km*1000 ;  "/,  m/sec 

vf_new  =  vf_ref  -  vf_atm;  "/,  m/sec 

change_v  =  (vf_new  -  v0)/v0*100;  "/,  percent  change  v 

'/,  Turn-off  the  control: 
u  =  [0  0]’; 

'/,  Call  ode45  to  propagate  the  satellite  EQM  forward  in  time  over  the  number 
'/,  of  periods  specified  for  the  post-optimization  analysis: 

[t_num2 , out2]  =  ode45 (<3cw2_eom, tspan,xO_new, options)  ; 

•/. -  PRINT  TO  SCREEN  - 

'/,  Output  the  total  control  used: 

disp(sprintf  ( ’  Integral  cost,  Total:  "/05 . 2f  ’  ,  costl)  )  ; 
disp(sprintf  ( ’Terminal  cost,  Total:  ”/02 . 7f  ’  ,  cost2)  )  ; 
disp(sprintf  (’Control  usage.  Sat  1:  ”/02 . 5f  m/sec ’  ,  control_totl)  )  ; 
disp(sprintf  (’Control  usage.  Sat  2:  ”/02 . 5f  m/sec ’,  control_tot2)  )  ; 
control_tot_init  =  control_totla  +  control_tot2a; 

disp(sprintf  (’Control  Check  (Total  from  orig  control):  '/,2.5f  m/sec’,... 
control_tot_init) ) ; 

disp(sprintf  ( ’Altitude  loss  (avg)  :  ”/,6.2f  m’ ,  altitude_loss2)  ) 
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disp(sprintf  ( ’Percent  Change  in  Mean  Motion:  "/,1 . 6f  ’  ,percent_change_n)  )  ; 
disp(sprintf ( ’Percent  Change  in  Relative  Velocity:  7,1 . 6f  ’  ,  change_v)  )  ; 
disp(sprintf ( ’Percent  Change  in  Atmospheric  Density:  7,1 . 6f  ’  ,  change_rho)  )  ; 

7,  Output  the  length  of  time  the  satellites  were  controlled  (specified  in 
7,  terms  of  orbital  periods,  but  convert  to  hours :  min:  sec)  : 
tp_hr  =  floor (t_fin/(3600) ) ;  diffl  =  t_f in-tp_hr*3600 ;  tp_min  = 
f loor (diff 1/60) ;  tp_sec  =  diff l-tp_min*60;  disp(’  ’); 

disp(sprintf  (’System  controlled  over  7,2. Of  hr,  7, 2. Of  min,  7,2. If  sec’,... 
tp_hr ,tp_min,tp_sec) ) ; 

7,  Output  the  terminal  constraints: 

disp(’  ’);  disp( ’Terminal  Constraints  from  NUMERICAL 
INTEGRATION: ’) 

disp(sprintf  ( ’Drift  Constraint:  7,0 . 5g’ ,  termc_drift_num)  ) 

disp(sprintf  ( ’Distance  Constraint:  7,4.  If  (m)  ’  ,termc_distance_num)  ) 

7. -  BEGIN  PLOTTING  - 

7,  Plot  the  results  of  the  optimization  (plot  of  controlled  portion)  : 

7,  Plot  the  control  history: 
figure (1)  subplot (2 , 1 , 1) 

plot (t/t_period,u_lin(l , : ) , ’k’ , t/t_period,u_lin(2 , : ) , ’k: ’ , ’LineWidth’ ,1.2) 
set(gca, ’XLim’ , [0,t_f inal] )  xlabel(’Time  (Orbital  Periods)’) 
ylabel ( ’Control  Input:  BC  (m~2/kg)’) 

7,  Convert  Figure  1  to  .  eps : 

filel  =  [foldername  ’ \control_input ’ ] ;  print (’ -fl ’,’ -depsc ’ ,filel) 

7,  Plot  the  radial/in-track  plane  for  the  numerical  integration: 
figure (2)  subplot (2 , 1 , 1) 

plot (x_num( : , 2) ,x_num( : , 1) , ’k’ ,x_num( : ,8) ,x_num( :,7),’k:’,0,0,’ko’,... 

0,0, ’kx’ , ’LineWidth’ ,1.2) 

x_min  =  min(min(x_num( : ,2) ) ,min(x_num( : ,8) ) ) ;  y_min  = 
min(min(x_num( : , 1) ) ,min(x_num( : ,7) ) ) ;  x_max  = 
max(max(x_num( : , 2) ) ,max(x_num( : ,8) ) ) ;  y_max  = 
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max(max(x_num( : , 1) ) ,max(x_num( : ,7) ) ) ; 

set (gca, ’ XLim’ , [1 . l*x_min, 1 . l*x_max] , ’ YLim’ , [1 . l*y_min, 1 . l*y_max] ) ; 
xlabel (’ In-track  (m) ’ )  ylabel ( ’Radial  (m) ’ )  legend (’ Satellite 
1’ , ’Satellite  2’) 

'/,  Convert  Figure  2  to  .  eps : 

filel  =  [foldername  ’ \f alling_f rame ’ ] ;  print ( ’ -f 2 ’ , ’ -depsc ’ ,f ilel) 

'/,  Plot  the  radial  and  in-track  time  histories  of  the  numerical  integration: 
figure(3)  subplot (2 , 1 , 1)  plot(t_num/t_period,x_num( : , 1) , ’k’ , . . . 

t_num/t_period,x_num( : ,7) , ’k: ’ , ’LineWidth’ ,1.2) 
set(gca, ’XLim’ , [0,t_f inal] )  xlabel(’Time  (Orbital  Periods)’) 
ylabel ( ’ Radial  (m) ’ ) 

subplot (2,1,2)  plot (t_num/t_period,x_num( :,2),’k’,... 

t_num/t_period,x_num( : ,8) , ’k: ’ , ’LineWidth’ ,1.2) 
set(gca, ’XLim’ , [0,t_f inal] )  xlabel(’Time  (Orbital  Periods)’) 
ylabel (’ In-track  (m) ’ ) 

'/,  Convert  Figure  3  to  .eps: 

filel  =  [foldername  ’\reconfig_time_history’] ; 
print ( ’ -f 3 ’ , ’ -depsc ’ , filel) 

figure (4)  subplot (2 , 1 , 1)  if  fig  ==  1, 

plot (out2 ( : , 2) , out2( : , 1) , ’k+ ’ ,out2( : ,8) , out2( : ,7) , ’k+’ , . . . 

0,0, ’k+> ,0,0, ’ko’ , ’LineWidth’ ,1.2) 
else,  plot ( out 2 ( : ,2) ,out2( : ,1) , ’k’ ,out2( : ,8) ,out2( :,7),’k:’,... 

0,0, ’k+> ,0,0, ’ko’ , ’LineWidth’ ,1.2) 
end  x_min  =  min(min(out2( : ,2) ) ,min(out2( : ,8) ) ) ;  x_max  = 
max(max(out2 ( : , 2) ) ,max(out2 ( : ,8) ) ) ;  y_min  = 
min(min(out2 ( : , 1) ) ,min(out2 ( : ,7) ) ) ;  y_max  = 
max(max(out2 ( : , 1) ) , max (out 2 ( : ,7) ) ) ; 

set (gca, ’ XLim’ , [1 . 2*x_min, 1 . 2*x_max] , ’YLim’ , [1 . 2*y_min, 1 . 2*y_max] ) ; 
xlabel (’ In-track  (m) ’ )  ylabel ( ’Radial  (m) ’ )  if  fig  ==  1 
axis  equal 

end 
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'/,  Convert  Figure  4  to  .  eps : 

filel  =  [foldername  ’ \post_inplane2 ’ ] ;  print (’ -f4’ ,  ’ -depsc ’ ,filel) 

'/,  Post-Optimization  time  history  plots: 
figure (5)  subplot (2 , 1 , 1)  if  fig  ==  1, 
plot  Ct_num2/t_period, out2 

t_num2/t_period, out2( : ,7) , ’k : ’ , ’ LineWidth’ ,1.2) 
else,  plot(t_num2/t_period,out2( : , 1) , ’k’ , . . . 

t_num2/t_period, out2( : ,7) , ’k: ’ , ’ LineWidth’ ,1.2) 
end  y_min  =  min(min(out2( : , 1) ) ,min(out2( : ,7) ) ) ;  y_max  = 
max(max(out2 ( : , 1) ) ,max(out2 ( : ,7) ) ) ;  if  fig  ==  1, 
set (gca, ’ YLim’ ,  [-0 .01,0.01]);  else , 

set(gca, ’YLim’ , [1 .4*y_min, 1 .4*y_max] ) ;  end  xlabel(’Time  (Orbital 
Periods)’)  ylabel ( ’Radial  (m) ’ ) 

subplot (2,1,2)  plot (t_num2/t_period, out2 ( : ,2) , ’k’ , . . . 

t_num2/t_period, out2( : ,8) , ’k: ’ , ’LineWidth’ ,1.2) 
y_min  =  min(min(out2( : , 2) ) ,min(out2 ( : , 8) ) ) ;  y_max  = 
max(max(out2 ( : , 2) ) ,max(out2 ( : ,8) ) ) ; 

set(gca, ’YLim’ , [1 .4*y_min, 1 .4*y_max] ) ;  xlabel(’Time  (Orbital 
Periods)’)  ylabel (’ In-track  (m) ’ ) 

'/,  Convert  Figure  5  to  .eps: 

filel  =  [foldername  ’ \post_time_hist ’] ; 

print ( ’ -f 5 ’ , ’ -depsc ’ , f ilel) 

figure (6)  subplot (2 , 1 , 1)  if  fig  ==  1 
yl_norm  =  norm(out2( : , 1)) ; 
y2_norm  =  norm(out2( : ,7)) ; 
plot (t_num2/t_period, out2( :,l),’k’,... 

t_num2/t_period, out2( : ,7) , ’k : ’ , ’ LineWidth’ ,1.2) 

else 

plot (t_num2/t_period, out2( : , 1) -d_sepf /4*sin(n0*t_num2) -out2(l , 1) , . . . 
’k’,  t_num2/t_period, . . . 

out2 ( : , 7)+d_sepf /4*sin(n0*t_num2)-out2(l , 7) , ’k : ’ , ’ LineWidth’ ,1.2) 
end  set (gca, ’YLim’ ,  [-0 . 008 ,0 . 008] ) ;  xlabel(’Time  (Orbital 
Periods)’)  ylabel ( ’Radial  Position  Deviation  (m) ’ ) 
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subplot (2 ,1,2)  if  fig  ==  1 


plot (t_num2/t_period, (out2 ( : , 2) -d_sepf /2) - . . . 

(out2( : ,8) +d_sepf /2) , ’k’ , ’LineWidth’ ,1.2) 
xlabelC’Time  (Orbital  Periods)’) 
ylabel( ’ In-Track  Deviation  (m) ’ ) 
elseif  fig  ==  2 

plot  (t_num2/t_period, (out2 ( : , 2) -d_sepf /2*cos (n0*t_num2) ) - . . . 

(out2( : ,8) +d_sepf /2*cos (n0*t_num2) ) , ’k’ , ’LineWidth’ ,1.2) 
xlabelC’Time  (Orbital  Periods)’) 
ylabel( ’ In-Track  Deviation  (m) ’ ) 

end 

'/,  Convert  Figure  6  to  .  eps : 

filel  =  [foldername  ’\post_hist_deviation’] ; 
print ( ’ -f 6 ’ , ’ -depsc ’ ,f ilel) 

'/,  If  we  decided  to  run  the  simulation  holding  the  relative  frame  fixed 
'/,  let’s  plot  the  results: 
if  fixed_frame  ==  1 

"/,  Plot  these  results: 
f igure(7) 
subplot (2,1,1) 

plot(x_new( : ,2) ,x_new( : , 1) , ’k’ .... 

x_new( : ,8) ,x_new( : ,7) , ’k: ’ , ’LineWidth’ ,1.2) 
x_min  =  min(min(x_new( : , 2) ) ,min(x_new( : ,8) ) ) ; 
y_min  =  min(min(x_new( : , 1) ) ,min(x_new( : ,7) ) ) ; 
x_max  =  max(max(x_new( : , 2) ) ,max(x_new( : ,8) ) ) ; 
y_max  =  max(max(x_new( : , 1) ) ,max(x_new( : ,7) ) ) ; 

set (gca, ’XLim’ , [1 . l*x_min, 1 . l*x_max] , ’YLim’ , [1 . l*y_min, 1 . l*y_max] ) 

xlabel( ’ In-track  (m) ’ ) 

ylabel( ’Radial  (m) ’ ) 

legend( ’ Satellite  1 ’,’ Satellite  2’) 

"/,  Convert  Figure  7  to  .eps: 

filel  =  [foldername  ’ \f ixed_frame ’ ] ; 

print ( ’ -f7’ , ’ -depsc ’ , filel) 
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%  Now,  we  should  propagate  this  final  state  forward  without  control  to 
"/,  show  that  it’s  equivalent  to  the  other, 
nm  =  length(x_new) ; 
xf_new  =  x_new(nm,:); 
u  =  [0 ;  0]  ; 

t_fin2  =  t_prop*t_period; 
tspan  =  [0  t_f  in2]  ; 

[tf_n,xf_n]  =  ode45(@cw2_eom, tspan, x0_new, options) ; 
else,  "/,  Do  nothing 
end 

'/,  Print  to  screen: 

disp(’  ’);  disp(’Run  completed.  Congratulations.’) 

toe  "/,  stop  clock 

return 

'/. -  NEW  FUNCTION:  CW2_E0M.M  - 

'/,  Function  cw2_eom.m 

'/,  Clohessy-Wiltshire  EOM  for  two  satellites  with  a  control  input. 

'/,  The  current  state  is  passed  as  input  to  the  function,  but  the  mean  motion 
'/,  (n)  ,  atmospheric  density  (rho)  ,  and  control  input  (u)  also  need  to  be 
'/,  passed  to  the  function. 

function  [dX]  =  cw2_eom(t,x)  global  nO  rho  u  vO 

'/,  Initialize  commonly  used  values: 
n  =  nO;  n2  =  n~2; 

'/,  Define  the  Earth  radius,  Gravitational  parameter,  and  J2: 

Re  =  6378.137;  7.  km 

mu  =  398600.4418;  ’/.  km~3/sec~2 

l  Define  the  A1  matrix  (6x6) : 


A1  =  [0,  0,  0,  1,  0,  0; 

0, 

0, 

0, 

0, 

1,  0;  0, 

0, 

0, 

t — 1 

O 

o 

3*n2,  0,  0,  0,  2*n, 

0; 

0, 

0, 

0, 

-2*n,  0, 

0; 

0, 

0 ,  -n2 , 
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'/,  Now  form  the  overall  A  matrix.  Since  the  mean  motion  corresponds  to  the 
'/,  reference  orbit,  A1  =  A2.  So  form  12x12  A  matrix  using  same  entries: 

A  =  [Al,  zeros(6,6);  zeros(6,6),  Al] ; 

'/,  Define  the  control  matrix  for  each  satellite  (6x1)  : 

B1  =  [0;  0;  0;  0;  -0 . 5*rho* (vO) *2 ;  0];  B2  =  [0;  0;  0;  0; 

-0.5*rho*(v0)~2;  0]; 

'/,  Now  form  the  overall  control  matrix  (12x2)  : 

B  =  [Bl,  zeros (6,1);  zeros (6,1),  B2] ; 

•/.  Define  the  EDM: 
dX  =  A*x  +  B*u;  return 

'/. -  NEW  FUNCTION:  HAM_E0M  - 

'/,  Function  ham_eom.m 

'/,  Subfunction  called  by  ’ode45.m’  to  propagate  the  Hamiltonian  E0M  forward 
'/,  in  time . 

function  [dX]  =  ham_eom(t,x)  global  nO  rho  u  vO 

'/,  Initialize  commonly  used  values: 
n  =  nO;  n2  =  n~2;  v06  =  v0“6;  rho2  =  rho"2; 

'/,  The  input  currently  contains  both  the  states  and  costates.  Split  the 
'/,  input  based  on  the  number  of  states  so  the  EQM  can  be  easily  defined. 
state_length  =  12;  xO  =  x(l : state_length) ;  lamO  = 
x(state_length+l :length(x)) ; 

'/,  Define  the  Al  matrix  (6x6)  for  the  state  of  1  satellite: 


Al  =  [0,  0,  0,  1,  0,  0; 

0, 

0, 

0, 

0, 

1,  0;  0, 

0, 

0, 

t — 1 

O 

o 

3*n2,  0,  0,  0,  2*n, 

0; 

0, 

0, 

0, 

-2*n,  0, 

0; 

0, 

0 ,  -n2 , 

'/,  Now  form  the  overall  A  matrix.  Since  the  mean  motion  corresponds  to  the 
'/,  reference  orbit,  Al  =  A2.  So  form  12x12  A  matrix  using  same  entries: 

A  =  [Al,  zeros(6,6);  zeros(6,6),  Al] ; 
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'/,  Define  the  control  matrix  for  each  satellite  (6x1)  : 

B1  =  [0;  0;  0;  0;  -0 . 5*rho* (vO) *2 ;  0];  B2  =  [0;  0;  0;  0; 

-0.5*rho*(v0)~2;  0]; 

'/,  Now  form  the  overall  control  matrix  (12x2)  : 

B  =  [Bl,  zeros (6,1);  zeros (6,1),  B2] ; 

'/,  Define  the  control  penalty  matrix  R: 

R  =  [0 . 25*rho2*v06 ,  0;  0,  0 . 25*rho2*v06] ; 

'/,  Define  the  Hamiltonian  dynamics: 

dx  =  A*x0  -  B*inv(R) *B’ *lam0 ;  dlambda  =  -A’*lam0; 

'/,  Form  the  output  vector: 
dX  =  [dx;  dlambda] ;  return 

'/. -  ATM0S_EXP.M  - 

'/,  B .  Haj  ovsky ,  30  Aug  06 

'/,  Exponential  Atmospheric  Model 

'/,  Accepts  as  INPUT  the  orbital  radius  of  the  satellite  in  km;  OUTPUTS  the 
'/,  atmospheric  density  for  that  altitude  based  on  the  mean  equatorial  radius 
'/,  of  Earth  (density  units:  kg/m~3)  . 

'/,  NOTE:  really  only  valid  for  altitudes  less  than  1000  km 
'/,  Taken  from  Vallado,  pg  537 
'/,  Distances  in  km 

function  [rho_exp] =atmos_exp(R) 

'/,  Term  Definitions: 

'/,  hO  =  Base  Altitude  (km) 

'/,  rhoO  =  Nominal  Density  (kg/m~3) 

'/,  H  =  Scale  Height  (km) 

'/,  Define  the  Earth  radius  to  calculate  the  altitude: 

Earth_rad  =  6378.137;  "/,  km;  mean  equatorial  radius 
r  =  R  -  Earth_rad; 
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if  r>=0  k  r<25,  hO  =  0;  rhoO  =  1.225;  H  =  7.249;  elseif  r>=25  k 
r<30,  hO  =  25;  rhoO  =  3.899e-2;  H  =  6.349;  elseif  r>=30  k  r<40,  hO 
=  30;  rhoO  =  1.774e-2;  H  =  6.682;  elseif  r>=40  k  r<50,  hO  =  40; 
rhoO  =  3.972e-3;  H  =  7.554;  elseif  r>=50  k  r<60,  hO  =  50;  rhoO  = 


1 . 057e-3 ; 

H  =  8.382; 

elseif 

r>=60 

k 

A 

~n] 

O 

hO  = 

60; 

rhoO  = 

3.206e-4; 

H  =  7.714; 

elseif 

Hi 

V 

II 

o 

k 

Hi 

A 

00 

o 

hO  = 

70; 

rhoO  = 

8 . 770e-5 ; 

H  =  6.549; 

elseif 

Hi 

V 

II 

00 

o 

k 

r<90 , 

hO  = 

80; 

rhoO  = 

1 . 905e-5 ; 

H  =  5.799; 

elseif 

r>=90 

k 

Hi 

A 

1-^ 

o 

o 

hO 

=  90 

rhoO 

3 . 396e-6 ; 

H  =  5.382; 

elseif 

r>=100 

k 

r<110 , 

hO  =  100; 

rhoO  = 

5 . 297e-7 ; 

H  =  5.877; 

elseif 

r>=110 

k 

r<120 , 

hO  =  110; 

rhoO  = 

9 . 661e-8 ; 

H  =  7.263; 

elseif 

r>=120 

k 

r<130 , 

hO  =  120; 

rhoO  = 

2.438e-8; 

H  =  9.473; 

elseif 

r>=130 

k 

r<140 , 

hO  =  130; 

rhoO  = 

8.484e-9; 

H  =  12.636; 

elseif 

r>=140 

k 

r<150, 

hO  =  140; 

rhoO  = 

3 . 845e-9 ; 

H  =  16.149; 

elseif 

r>=150 

k 

r<180, 

hO  =  150; 

rhoO  = 

2 . 070e-9 ; 

H  =  22.523; 

elseif 

r>=180 

k 

r<200, 

hO  =  180; 

rhoO  = 

5.464e-10; 

H 

= 

29.740; 

elseif 

r>=200 

k 

r<250 , 

hO  = 

200; 

rhoO  = 

2 . 789e-10; 

H 

= 

37.105; 

elseif 

r>=250 

k 

r<300 , 

hO  = 

250; 

rhoO  = 

7 . 248e-ll ; 

H 

= 

45.546; 

elseif 

r>=300 

k 

r<350 , 

hO  = 

300; 

rhoO  = 

2.418e-ll; 

H 

= 

53.628; 

elseif 

r>=350 

k 

r<400 , 

hO  = 

350; 

rhoO  = 

9 . 518e-12 ; 

H 

= 

53.298; 

elseif 

r>=400 

k 

r<450 , 

hO  = 

400; 

rhoO  = 

3 . 725e-12 ; 

H 

= 

58.515; 

elseif 

r>=450 

k 

r<500 , 

hO  = 

450; 

rhoO  = 

1 . 585e-12 ; 

H 

= 

60.828; 

elseif 

r>=500 

k 

r<600 , 

hO  = 

500; 

rhoO  = 

6 . 967e-13; 

H 

= 

63.822; 

elseif 

r>=600 

k 

r<700 , 

hO  = 

600; 

rhoO  = 

1.454e-13; 

H 

= 

71.835; 

elseif 

r>=700 

k 

Hi 

A 

00 

o 

o 

hO  = 

700; 

rhoO  = 

3.614e-14; 

H 

= 

88.667; 

elseif 

r>=800 

k 

Hi 

A 

CO 

O 

O 

hO  = 

800; 

rhoO  = 

1 . 170e-14; 

H 

124.64; 

elseif 

r>=900 

k 

O 

O 

O 

tH 

V 

hO 

=  900 

rhoO  = 

5 . 245e-15 ; 

H 

= 

181.05; 

elseif 

r>=1000. 

hO  =  1000; 

rhoO 

= 

3 . 019e-15 ; 

H 

= 

268.00; 

else , 

I1RR0R  (  ’  Che  ck  Atmo  sphe 

ric  Density 

file!!  Correct  radius  not  found!’)  end 
'/,  Output  the  density  measurement : 
rho_exp  =  rhoO*exp(-(r-hO)/H) ;  return 

'/. -  END  SCRIPT  - 
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